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ABSTRACT 

This  thesis  is  an  extension  of  the  theory  of 
discrete  scalar  time  series  analysis  to  multivariable 
processes.  This  extension  is  facilitated  by  expanding 
the  algebra  of  polynomial  matrices  (matrices  with  poly¬ 
nomial  elements ) . 

Multivariable  processes  may  have  a  multiplicity 
of  eitb  r»  the  independent  or  the  dependent  variable.  Such 
processes  are  called  multi-dimensional  or  multi-channel, 
respectively.  All  multi-dimensional  processes  may  be 
formally  mapped  into  matrix  notation.  Once  this  mapping 
is  made  the  properties  of  all  multivariable  linear  opera¬ 
tors  and  autocorrelations  can  be  studied  in  terms  of  the 
polynomial  matrices  that  represent  their  z-transf orms . 

Polynomial  matrices  can  be  decomposed  into  three 
related  forms:  the  spectral  factorization,  the  Smith- 
McMillan  canonical  form,  or  the  Robinson  canonical  form. 
Each  of  these  representations  leads  to  the  concept  of  an 
invertible  or  minimum  delay  wavelet. 

The  algorithms  for  finding  the  spectral  factori¬ 
zation  and  for  finding  the  Smith-McMillan  canonical  form 
can  be  extended  to  provide  an  analytic  factorization  of  a 
multi-channel  autocorrelation  in  term  of  invertible  wave¬ 
lets.  In  addition  the  autocorrelation  may  be  approximately 
factored  by  a  recursive  least-squares  algorithm,  or  by  a 
projection  technique. 

Of  the  factorization  methods  available,  the  re¬ 
cursive  algorithm  is  the  most  efficient  and  is  therefore 
extended  to  include  the  more  general  problem  of  signal 
shaping  in  the  presence  of  noise. 

Finally,  as  an  illustration,  the  problem  of 
designing  a  finite  optimum  two-dimensional  band-pass,  band- 
reject  filter  is  solved  and  the  characteristics  of  a  few 
particular  realizations  of  such  filters  are  presented. 

Thesis  Supervisor:  Stephen  M.  Simpson,  Jr. 

Title:  Lecturer  in  Geology  and  Geophysics 
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1.  INTRODUCTION 


Geophysics  may  be  viewed  as  the  study  of  the 
properties  of  the  earth  by  the  interpretation  of  signals 
that  are  affected  by  the  structure  of  the  earth.  These 
signals  may  be  of  almost  any  conceivable  type  --  seismic, 
tidal,  electric  current,  electromagnetic,  or  light  —  and 
may  have  a  wide  variety  of  sources.  In  each  case  the  com¬ 
plexity  of  the  media  that  modulates  the  signals  will  intro 
duce  ncise  into  the  system  (we  define  noise  as  any  portion 
of  a  signal  which  does  not  contain  information  that  we  de¬ 
sire).  In  addition,  the  information  may  be  difficult  to 
interpret  because  the  signal  shapes  are  difficult  to 
recognize. 


The  idea  of  applying  the  concepts  of  statistical 
analysis  to  signal  interpretation  has  become  widely  accepted 
during  the  last  decade.  A  large  portion  of  this  analysis 
has  taken  the  form  of  applying  linear  filters  to  incoming 
data  to  enhance  and  shape  the  desired  information.  Because 
of  its  versatility,  the  least-squares  optimum  (Wiener)  fil¬ 
ter  was  frequently  applied.  However,  a  problem  arose  in 
the  computation  of  such  filters  for  geophysical  applica¬ 
tions.  Geophysical  signals  are  usually  multivariable;  that 
is,  the  signals  are  characterized  by  having  more  than  one 
independent  variable  (dimension)  or  by  having  more  than  one 
dependent  variable  (channel).  In  most  cases,  to  adequately 
process  such  multivariable  signals  one  should  use  multi- 
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variable  filters.  However,  the  solution  of  the  discrete 
least-squares  filter  problem  is  a  set  of  simultaneous  equa 
tions,  with  one  equation  for  each  coefficient.  Thus  the 
magnitude  of  the  problem  quickly  overloaded  the  capacity 
of  even  the  largest  computers  available.  This  limitation 
on  the  size  of  possible  filters  greatly  restricted  their 
usefulness. 


The  computational  problem  was  reduced  by  an  order 
of  magnitude  in  storage  space  and  execution  time  by 
Robinson  (1963a)  when  he  was  able  to  extend  a  recursive 
method  introduced  by  Levinson  (19^7)  to  multivariable 
filter  generation.  This  development  has  led  to  greatly 
renewed  Interest  in  the  applications  of  optimum  filtering. 

The  crucial  step  in  optimum  least-squares  filter 
design  is  the  factorization  of  the  autocorrelation  of  a 
process.  This  factorization  is  the  problem  to  which  this 
thesis  is  addressed. 

There  are  now  four  known  techniques  for  factoring 
the  autocorrelations  of  multivariable  processes.  As  indi¬ 
cated  above,  the  least-squares  approximate  factorization 
has  been  known  for  some  time  (for  example,  see  Wadsworth, 
et  al.,  1963).  However,  before  the  discovery  of  the  re¬ 
cursive  computation  algorithm,  it  was  not  considered  to  be 
useful.  In  fact  it  was  this  consideration  that  led  Wiener 
and  Masanl  (1957  and  1956)  to  develop  a  projection  technique 
of  approximate  factorization.  Experience  now  shows  that 


12 


this  technique  is  not  competitive  with  the  recursive  method. 
Quenoullle  (1957)  presented  an  analytic  factorization  al¬ 
gorithm  which,  when  placed  upon  a  rigorous  mathematical 
basis,  has  proven  to  be  a  very  valuable  tool  for  under¬ 
standing  and  manipulating  multivariable  time  series  and 
autocorrelations.  Another  analytic  factorization  method 
is  developed  here  based  upon  the  Smith  canonical  form  for 
polynomial  matrices  following  a  similar  development  by 
You la  (1961) .  Neither  of  these  analytic  methods  are  compu¬ 
tationally  competitive  with  the  least-squares  recursive 
algorithm  although  they  are  invaluable  for  instilling 
theoretical  insight  into  the  factorization  problem. 

All  of  the  factorization  schemes  that  are  consid¬ 
ered  are  stated  for  discrete  processes  with  finite  auto¬ 
correlations.  Since  these  factors  (which  we  call  wavelets) 
are  also  finite  they  are  members  of  the  Hardy  class 
(Wiener  and  Masanl,  1957»  P.  113) •  Because  we  are  dealing 
with  finite  wavelets  we  are  able  to  obtain  specific  results 
which  are  of  a  more  constructive  nature  than  those  found 
in  some  more  generalized  approaches.  As  such  this  thesis 
may  be  considered  as  a  complement  to  recent  works  on  stoch¬ 
astic  processes  such  as  He Ison  and  Lowdenslager  (1958), 
Robinson  (1962),  and  Wiener  and  Masanl  (1957  and  1958). 

This  thesis  tnen  is  primarily  an  examination  and 
evaluation  of  the  methods  of  factoring  multivariable  auto¬ 
correlations.  Prom  another  point  of  view,  however,  it  may 
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be  thought  of  as  a  treatise  on  polynomial  and  rational 
matrices,  that  is,  on  matrices  whose  elements  are  either 
polynomial  or  rational.  This  is  a  subject  that  has  re¬ 
ceived  surprisingly  little  attention  in  the  literature. 

For  this  reason  it  is  given  a  rather  thorough  development 
here  in  the  first  three  chapters. 

The  final  chapters  are  devoted  to  an  expansion 
of  the  least-squares  approximate  factorization  to  the 
calculation  of  filters  with  specified  noise  suppressing 
and  signal  shaping  properties.  Computational  examples  are 
included  that  illustrate  some  of  the  forms  that  such 
computations  may  take. 

The  presentation  that  follows  assumes  a  basic 
knowledge  of  scalar,  i.e.  single-variable,  time  series 
analysis  (see  Lee,  I960;  Robinson,  1962;  Wiener,  1949;  or 
Whittle,  1963).  Most  of  the  primary  ideas,  such  as  wave¬ 
lets,  all-pass  systems,  minimum  phase,  minimum  delay, 
convolution,  autocorrelation,  and  predictive  decomposition 
are  reviewed  briefly  when  they  are  first  encountered  but 
are  not  developed  rigorously.  The  material  here  is  not 
Intended  to  be  a  review  of  time  series  analysis,  but  is 
intended  to  be  an  extension  of  the  concepts  of  scalar  time 
series  to  multivariable  time  series.  On  the  other  hand, 
much  of  the  detail  considered  is  not  necessary  for  an  over¬ 
all  grasp  of  multivariable  time  series  analysis.  Thus  the 
reader  who  is  unfamiliar  with  the  subject  may  profitably 
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skip  over  several  sections.  These  sections  include  4.221 
(the  details  of  the  factorization  of  an  elementary  auto¬ 
correlation  matrix),  4.223  (the  Smi th-McMillan  factorization 
technique),  4.232  (the  Wiener-Masani  approximate  factori¬ 
zation  by  projections),  and  5.12  (the  details  of  the  re¬ 
cursive  algorithm  for  least-squares  filters). 
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2.  DEFINITIONS  AND  NOTATION 


Processes  may  have  multiple  independent  or  depen¬ 
dent  variables.  In  sections  3.4  and  4.3  a  technique  is 
developed  for  mapping  processes  with  several  independent 
variables  into  a  form  with  several  dependent  variables  and 
only  one  independent  variable.  This  mapping  is  given  in 
order  to  simplify  the  analysis  and  notation  of  the  factori¬ 
zation  problem.  However,  there  are  important  differences 
between  these  two  representations  that  should  be  recalled 
when  applications  are  made  of  factorization.  This  chapter 
is  devoted  primarily  to  an  examination  of  these  differences. 
In  addition,  a  few  general  notational  questions  are  exam¬ 
ined. 


2.1  Dimensionality  of  Processes 

A  dimension  is  defined  as  a  measurable  extent. 

In  this  thesis  the  number  of  dimensions  of  a  process  will 
indicate  the  number  of  orthogonal  measurable  directions, 
i.e.,  the  number  of  Independent  variables.  Most  processes 
that  have  been  considered  in  communication  theory  and  in 
economic  analysis  are  one-dimensional  time  series.  However, 
in  geophysics  higher  dimensioned  processes  are  often  en¬ 
countered  that  may  or  may  not  have  a  time-like  dimension. 
For  example,  the  output  of  a  single  vertical  seismometer 
is  a  one-dimensional  time  process.  The  output  of  a  linear 
row  of  seismomenters  is  a  two-dimensional  process  —  one 
time  dimension  and  one  space  dimension.  On  the  other  hand 
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(neglecting  small,  higher  order  effects)  the  acceleration 
of  gravity  at  each  of  these  seismometer  locations  repre¬ 
sents  a  one-dimensional  spatial  process. 

In  nearly  '\11  of  our  analysis  we  shall  assume 
that  one  of  the  dimensions,  or  directions,  of  a  process  is 
a  preferred  (time-like)  direction.  We  do  this  for  several 
reasons,  (l)  In  many  processes  there  actually  exists  a 
preferred  direction.  It  is  only  natural  to  take  advantage 
of  the  physical  significance  of  this  direction.  (2)  The 
use  of  vector  notation  greatly  simplifies  the  representation 
of  processes  with  a  preferred  direction.  (3)  Present  digi¬ 
tal  computers  have  one-dimensional  storage  memories.  Thus 
when  a  process  is  mapped  into  a  computational  scheme,  we 
must  necessarily  choose  a  preferred  direction. 

Whittle  (195**  >  P*  **3*0  has  pointed  out  that  there 
is  a  basic  difference  between  a  preferred  direction  that 
has  time,  or  time-like,  physical  significance  and  a  direc¬ 
tion  that  is  chosen  merely  for  notational  purposes.  A  time¬ 
like  direction  is  inherently  one-sided.  That  is,  the 
state  of  a  process  at  any  time  can  be  dependent  only  upon 
past  values  of  the  process.  However,  purely  spatial  proc¬ 
esses  are  usually  not  one-sided.  This  distinction  is  im¬ 
portant  when  designing  operators  for  processes. 

The  importance  of  the  preferred  direction  is  em¬ 
phasized  when  we  define  the  geometrical  structure  of  the 
sampling  of  the  independent  variables  of  a  process.  We  can 
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think  of  this  structure  as  an  array  of  sample  points  in  a 

multi-dimensional  space.  For  nearly  all  applications  we 

will  restrict  these  points  to  being  equally  spaced  along 

straight  lines  that  are  parallel  to  the  preferred  direction. 

This  is  equivalent  to  saying  that  the  process  will  have 

equal  digitization  ‘ncrements  in  the  preferred  direction 

% 

and  fixed  sampling  Ltions  in  the  other  directions.  Al¬ 
though  the  digitization  increment  is  fixed,  the  sampling 
instants  for  the  various  positions  need  not  be  in  phase. 

For  most  applications  in  this  paper,  we  will  re¬ 
quire  that  the  lines  form  regular  patterns  in  the  other 
dimensions.  The  simplest,  and  most  useful,  pattern  is 
that  of  a  rectangular  grid.  However,  other  patterns 
(triangles,  parallelograms,  hexagons,  and  combinations  of 
these  in  higher  dimensions)  are  frequently  encountered. 

2.2  Order  of  Processes 

Processes  may  have  multiple  dependent  variables 
as  well  as  multiple  independent  variables.  In  general, 
the  dependent  variables  need  not  have  any  dimensional 
relationship.  For  example,  one  variable  may  represent  the 
electric  field  while  another  may  represent  the  magnetic 
field.  The  order  of  a  process  is  the  number  of  dependent 
variables  that  represent  a  process  at  each  point  in  space. 
Thus,  a  linear  array  of  3  component  seismometers  would  be 
a  3rd  order,  2-dimensional  process. 


lh 


Throughout  this  paper  we  will  refer  to  processes 
of  order  greater  than  1  as  multi-channel  or  matrix-valued. 
The  latter  designation  stems  from  the  fact  that  we  will 
use  a  matrix  representation  to  group  the  variables  of  a 
process . 


The  one-dimensional,  multi-channel  process  is  of 
special  interest  since  its  configuration  best  reflects  the 
Importance  of  the  preferred  direction.  This  fact  sometimes 
prompts  us  to  view  each  of  the  space  samples  of  a  multi¬ 
dimensional  process  as  one  channel  for  a  higher-ordered 
multi-channel  process.  Thus,  a  linear  array  of  12  three- 
component  seismometers  might  be  viewed  as  a  3rd  order,  2- 
dlmensional  process,  or,  viewing  each  seismometer  as  pro¬ 
viding  a  separate  time  series  (channel),  we  may  view  this 
as  a  36th  order,  one-dimensional  process. 

Even  though  a  mapping  from  a  multi-dimensional 
process  to  a  higher-ordered,  one-dimensional  process  is 
possible,  the  basic  differences  between  these  representa¬ 
tions  should  be  emphasized.  First,  we  usually  think  of  a 
discrete,  multi-dimensional  process  as  a  manifestation  of 
a  continuous  function.  Thus,  it  is  possible  to  approximate 
values  between  the  digitization  positions  by  some  form  of 
Interpolation.  Second,  in  a  multi-dimensional  process  we 
can  think  of  extending  the  space  dimensions  to  infinity. 

The  formal  structure  of  a  multi-channel  process  allows 
neither  of  these  possibilities. 
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2 . 3  Subscript  Notation 

Subscripts  will  be  used  to  indicate  the  variables 
of  a  process.  In  general,  there  will  be  two  groups  of  sub¬ 
scripts.  The  first  group  will  refer  to  Indexing  of  the 
independent  variables;  the  second  group  will  refer  to  in¬ 
dexing  of  the  dependent  variables.  We  will  adopt  the  con¬ 
vention  that  the  first  subscript  in  the  first  group  will 
always  stand  for  the  preferred  direction.  Thus  a  component 
of  a  process  X  may  be  referred  to  as 

i  1 

x,  . .  . , 

or,  if  the  preferred  direction  refers  specifically  to  time, 
it  will  be  written 

(^t  i  i  )  b-  ic 

u, x. ,  . . . ,  x^ 

Since  matrices  are  at  most  2-dimensional,  the  sec''  d  group 
will  have  at  most  2  indices.  We  will  always  consider  that 
the  first  of  these  2  indices  will  be  the  row  index,  and 
that  the  second  will  be  the  column  index. 

In  order  to  simplify  our  writing  we  shall  adopt 
a  vector  notation  for  the  subscripts 

1  =  •••»  1 n ) 

k  =  (k-^kg) 

so  that  the  process  may  also  be  referred  to  as 

(*i>k 


Frequently  It  is  desirable  to  order  the  spatial 
sampling  positions  (i.e.  the  sampling  positions  in  the 
non  preferred  direction)  sequentially.  Thus  we  may  use  one 
subscript  for  all  spatial  variables: 

This  subscript  takes  on  a  different  value  for  each  sampling 
position.  Finally,  for  much  of  our  work  we  will  be  con¬ 
cerned  only  with  the  dimensional  indices  and  will  suppress 
the  matrix  indices  and  the  parentheses. 


2.4  Flow  Diagram  Notation 

The  important  decomposition  and  factorization 
theorems  in  the  following  chapters  are  illustrated  by  flow 
diagrams.  In  general,  these  diagrams  are  self  explanatory, 
however,  a  description  of  some  of  the  conventions  used  will 
facilitate  their  interpretation. 

1.  Square  boxes  indicate  operations. 

2 .  Rounded  boxes  illustrate  results  of 

operations . 

3.  Sol'd  lines  between  boxes  indicate  the 

primary  lines  of  logical  flow  as  well 

as  transference  of  data  between  steps. 

4.  Dotted  lines  between  boxes  indicate  only 

the  transference  of  data  between  steps. 

5.  Boxes  drawn  with  heavy  lines  indicate  the 

beginning  and  the  ending  of  the  algorithm. 
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3.  STRUCTURE  OF  DISCRETE  LINEAR  OPERATORS 


The  operators  that  we  consider  are  finite  moving 
average  devices  that  may  be  represented  by  the  diagram 


Input 


Linear 

Operator 

Output 


If  the  input  is  a  spike  (a  delta  function  appropriate  to 
the  Geometry  involved)  then  the  output  is  a  wavelet  with 
real  coefficients  which  completely  describes  the  properties 
of  the  linear  operator.  In  fact  the  output  of  the  linear 
operator  for  a  general  input  is  just  the  convolution  of 
the  input  with  the  wavelet. 


In  this  chapter  we  will  study  those  characteris¬ 
tic  properties  of  a  wavelet  by  which  it  may  be  classified. 
The  approach  used  here  is  to  factor  a  wavelet  into  simpler 
components  and  then  to  use  the  properties  of  these  compo¬ 
nents  to  delineate  the  classification  of  the  wavelet.  The 
complexity  or  existence  of  the  factorization  is  the  key 
problem.  In  the  scalar  one-dimensional  case  there  is  a 
unique  natural  factorization  from  which  the  general  proper' 
ties  are  easily  deducible.  In  the  matrix-valued  one¬ 
dimensional  case  there  are  a  multiplicity  of  such  factori¬ 
zations.  In  the  multi-dimensional  case  there  is  no 
natural  factorization.  Thus,  our  treatment  for  these 
cases  will  vary  markedly. 


3,1  z- Transform 


The  z-transform  of  a  discrete  finite  wavelet  is 


defined  simply  as  the  quasipolynomial 


1 

a(z)  =  £  a .  z  ,  -  <*><  N  <  M  <  °° 

i=N  1 


whose  coefficients  a^  are  the  values  of  the  wavelet  at 
the  i~  sample  time.  For  the  general  multi-dimensional 
processes  we  have 


a(z,zx» 


M. 


zr)  = 


^  (ai  i  ^ 

i .=N  .  1,11'  ***'  1n 


)z 


i  i- 


■J""  J 


-  oo  <  Nj  <  Mj  <  oo  . 

A  quasipolynomial  a(z)  may  always  be  transformed  into  a 
polynomial  by  multiplying  it  by  the  proper  pc./er  of  z  . 

The  z-transform  of  a  wavelet  will  be  indicated  specifically 
by  writing  the  wavelet  as  a  function  of  z  as  indicated 
above . 


Two  important  properties  of  the  z-transform  will 
be  exploited  frequently: 

1.  Convolution  in  the  time-space  domain  corres¬ 
ponds  to  multiplication  in  the  z  domain. 

2.  The  z-transform  evaluated  on  the  unit  circle, 
z  =  e-^11’,  corresponds  to  the  Fourier  trans¬ 


form  of  the  wavelet 


Much  of  the  analysis  in  this  chapter  is  based 
on  the  algebra  of  quasipolynomials  that  corresponds  to 
the  z-transforms  of  wavelets. 

3.2  One-Dimensional  Scalar  Wavelets 

One-dimensional  scalar  wavelets  of  the  Hardy 
class  have  been  treated  extensively  in  the  literature, 

(Wold,  1938;  Wiener  and  Masanl ,  1957;  Robinson,  1962; 
Whittle,  1963;  Robinson  and  Treltel,  196^)  and,  therefore, 
the  treatment  here  will  be  brief  and  heuristic. 

3.21  Spectral  Decomposition 

Let  us  consider  the  one-sided  wavelet 
a0'  alf  an 

The  z-transform  of  this  wavelet 

a(z)  =  ao  +  al  z  +  • • •  +  ap  2 

can  be  factored,  according  to  the  fundamental  theorem  of 
algebra,  into  the  form 

a(z)  =  Uq(1  -  a j  z)  ...  (1  —  z) 

where  1  /a  ^  i  =  1,  ...,  n  are  the  zeros  of  the  polynomial 
a(z)  .  These  roots,  1  /a^,  are  generally  complex  but  since 
the  coefficients  of  a(z)  are  real,  the  roots  must  occur 
in  complex  conjugate  pairs. 


c. 


L 


3.22  Invertibility 


Definition  3.2-1.  A  one-sided  wavelet  a(z)  is  Slid  to  be 
invertible  if  there  exists  a  one-sided  wavelet  a  ^(z) 
such  that  a(z)  a  ^(z)  =  1  . 

The  condition  that  the  Taylor  expansion  of 
l/a(z)  will  converge  is  that  a(z)  has  no  zeros  inside 
the  unit  circle.  Thus  if  Jl/a^J  >1  i  =  1,  . ..,  n  then 
a(z)  is  invertible. 

Jury  (1964)  reviews  several  simple  techniques  for 
testing  for  the  invertibility  of  a  wavelet.  One  of  the 
simpler  conditions  involves  polynomial  divisions  to  find 
the  number  of  roots  inside  the  unit  circle.  The  procedure 
begins  by  performing  the  division 

a (z)  _  a1(z) 

zn  a(l/z)  0  zn  a(l/z) 

where  a^Cz)  is  the  remainder.  Then  we  find  the  other 
qi  i  =  1,  . n  -  2  according  to 

ax(z)  ai+l(z) 

9i  =  z"-1  H(l/z)  '  z^1  a^l/z) 

Now  the  number  of  roots  inside  the  unit  circle  is  equal  to 
the  number  of  products  which  are  negative,  where 

is  defined  as 
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3.23  Robinson  Canonical  Form  o.nd  All-Pass  Systems 


Theorem  3.2-1  (Robinson  Canonical  Form) .  Any  wavelet 

a(z)  can  be  uniquely  represented  in  the  Robinson  Canonical 

Form 

a(z)  =  p(z)  a0(z) 

where  aQ(z)  is  invertible  and  p(z)  is  an  all-pass 
system. 


Let  us  review  a  few  properties  of  all-pass 

systems . 


Theorem  3.2 -2 .  An  all-pass  discrete  system  has  unit  gain 
at  all  frequencies,  i.e.  |p(e”i"l)|  =  1  for  all  real  uj 


Theorem  3.2-3.  An  all-pass  system  is  trivial  if 
for  all  zj  that  is,  if  p(z)  is  constant. 


1 


Theorem  3.2-4 .  The  inverse  system  to  a  non-trivial  all¬ 
pass  system  is  not  one-sided. 


The  invertible  factor  aQ(z)  is  completely  det 
ermined  by  the  amplitude  spectrum  of  a(z)  (see  4.122). 

3.24  Delay 


The  delay  of  a  one-sided  wavelet  a^  is  a  measure 
of  how  the  operator  redistributes  the  energy  of  an  input 
process  in  forming  the  output.  It  may  be  defined  in  terms 
of  the  partial  energy 


dC 


Robinson  (1962)  has  proven  the  following  Minimum 


Delay  Theorem. 

Theorem  3,2-5  (Minimum  Delay) .  The  delays  of  the  set  of 
wavelets  a^(z)  which  have  the  same  invertible  Robinson 
canonical  form  aQ(z)  are  greater  than  or  equal  to  the 
delay  of  aQ(z)  .  Equality  holds  if  and  only  if  the  all¬ 
pass  system  p(z)  is  trivial.  That  is,  the  partial  ener¬ 
gies  obey  the  relation 

k  p  k  p 

T.  af  ,  <  T.  aZ  .  for  all  k 

j=0  1,J  j=0  U'J 

where  i  is  the  wavelet  index  and  j  is  the  time  index. 


3.25  Phase 


The  Fourier  transform  of  a  wavelet  yields  fre¬ 
quency  information  about  the  outputs  of  an  operator  with 
respect  to  the  inputs.  This  information  is  presented  in 
the  form  of  an  amplitude  change  and  a  phase  lag. 


If  we  examine  the  Fourier  transform  of  a  wavelet 


/  —  itn N  ,  -id)  ,  -i2i»  ,  ,  -inui 

a (e  )  =  an  +  an  e  +  a2  e  +  . . .  +  an  e 


=  a(«>)  e 


1 

iCD  (tU  ) 


we  see  ^hat  the  polar  representation  leads  to  the  concept 
of  a  phase  lag  characteristic  -m(t»)  . 


We  a  re  now  in  a  position  to  formulate  the 


minimum-phase  theorem: 

Theorem  3*2-6  (Minimum  phase) .  The  phase-lags  of  the  set 
of  wavelets  a^(z)  which  have  the  invertible  Robinson 
canonical  form  a^(z)  are  greater  than  or  equal  to  the 
phase-lag  of  aQ(z)  .  Equality  holds  if  and  only  if  the 
all-pass  system  p(z)  is  trivial.  Furthermore,  the 
phase-lag  difference  is 

^(0)  -  ^(tt)  =  m^rr 

where  m^  is  the  number  of  zeros  of  a^(z)  that  are  in¬ 
side  the  unit  circle. 

An  interesting  result  that  follows  directly  from 
the  theorem  above  is 

Corollary  3.2-6 .  The  cosine  transform  a (cos  o)  of  a 
wavelet  is  non-negative  if  the  wavelet  a(z)  has  no  zeros 
inside  the  unit  circle  (i.e.  is  minimum  phase  or  minimum 
delay)  and  if  the  wavelet  is  normalized  so  that  a(l)  >  1  „ 
The  number  of  zeros  of  a(z)  inside  the  unit  circle  is 
equal  to  the  number  of  zero  crossings  of  the  cosine- 
transform  a  (cos  u> )  , 

The  proofs  to  both  the  Minimum  Phase  Theorem  and 
its  Corollary  follow  directly  from  examining  the  nature  of 
the  definition  of  phase  (Robinson  and  Treltel,  1964). 

Figure  3.2-1  illustrates  the  behavior  of  the 


l,j  11  y  *  ■ 


phase-lag  curve  for  a  3-term  wavelet  for  various  positions 
of  the  zeros  near  the  unit  circle.  In  this  case  the  zeros 
were  placed  on  the  imaginary  axis  so  that  the  discontinuity 
for  the  middle  curve  lies  at  ««  =  rt/2  .  Because  of  this 
discontinuity  we  may  interpret  the  wavelet  either  as  mini¬ 
mum  phase  or  maximum  phase. 

3.3  One -Dimensional  Matrix-Valued  Wavelets 

Various  aspects  of  matrix-valued  wavelets,  or 
polynomial  matrices,  have  been  treated  by  a  number  of 
authors.  This  section  will  review  in  some  detail  many  of 
their  important  results  as  well  as  extend  the  theory  in 
certain  areas. 


3.31  Polynomial  Matrix  Notation 


Let  us  begin  by  reviewing  the  basic  notation  and 
terminology  used  in  describing  polynomial  matrices. 


Let  A  be  an  arbitrary  matrix.  Then: 


A' 


denotes  transpose 


denotes  complex  conjugate 
denotes  complex  conjugate  transpose 
denotes  inverse 


Det  A 

or  j  A  j  denotes  determinant  of  A 

\dj  A  denotes  adjugate  of  A.  (The  ad Jugate  of  A 
is  the  transposed  matrix  of  cofactors  of  A. 
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Note  that  Adj  A /bet  A  =  A-1  if 

UtA/O  .) 

A  diagonal  matrix  A  with  diagonal  terms  a,,  a„ ,  . ..,  a 

l  u  n 

is  written  as  A  =  dlag  £a^,  a-,,  ...»  arJ  •  Column  vectors 
are  represented  by  x>  1L>  etc.,  or  in  the  alternative 
fashion  x  =  (x^,  x^,  . ..,  x  )•  whenever  it  is  desirable 
to  indicate  the  components  explicitly.  The  symbols  1  or  I, 
0  .  and  0  represent  the  n  x  n  Identity  matrix,  the 
n-component  zero  vector  and  the  n  x  m  zero  matrix. 

A  matrix  A(z)  is  polynomial  or  quasipolynomial 
if  each  of  its  elements  is  a  polynomial  or  quasipolynomial 
in  z  .  A(z)  is  rational  if  each  of  its  elements  is  the 
ratio  of  two  polynomials  or  quasipolynomials  in  z  . 

A (z )  is  said  to  be  real  if  A(z)  =  A(z)  . 

Unless  stated  otherwise,  all  matrices  considered  here  will 
be  real. 


The  non-negative  integer  r(A)  is  the  rank  of 
the  rational  matrix  A(z)  for  a  given  value  of  z  if 

(1)  there  exists  at  least  one  subminor  of  order 
r  which  does  not  vanish  identically,  and 

(2)  all  minors  of  order  >  r  vanish  identi¬ 
cally. 

The  rank  of  an  n  x  n  matrix  A(z)  is  the  same  for  all 
z  except  for  a  finite  set  of  points  z^,  1=1,  ...,  p 

in  the  z  plane  at  which  the  rank  may  decrease.  These 
points  are  known  as  the  latent  zeros  of  the  matrix  A(z) 
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(see  section  3*32;.  The  maximum  number  of  latent  zeros  for 
an  n  x  n  matrix  A(z)  is  p  =  n  m  where  m  is  the  maxi¬ 
mum  number  of  zeros  in  any  quasipolynomial  element  of 
A(z)  .  .If  p  <  n  m  the  matrix  A(z)  is  called  degenerate . 

A  nonsquare  matrix  does  not  possess  an  inverse  in 

the  ordinary  sense.  However,  it  may  have  either  a  right 

or  left  Inverse .  Thus,  if  A  is  m  x  n,  A  possesses  a 

right  inverse  A-”1,  su  h  that  A  A-^  =1  if  and  only  if 

m 

m  <  n  and  r (A )  =  m  . 

An  elementary  polynomial  matrix  is  a  polynomial 
matrix  possessing  either  a  right  or  left  polynomial  inverse. 
A  square  matrix  A(z)  is  elementary  if  and  only  if  its 
determinant  is  independent  of  z  and  non-zero. 

A (z )  is  analytic  in  a  region  of  the  z  plane 
if  all  of  its  elements  are  analytic  in  this  region. 

The  point  zQ  is  a  pole  of  A(z)  if  son  ele¬ 
ment  of  A(z)  has  a  pole  at  z  =  zQ  . 

If  zQ  is  a  pole  of  the  rational  matrix  A(z), 
each  element  of  A  may  be  expanded  in  partial  fractions 
and  after  collecting  all  those  terms  having  poj.es  at  zQ 
there  is  obtained  for  zQ  /  co 

A(z)  =  (z  -  Zq )  Ak  +  (z  -  zQ)  Ak-1  +  ...  + 

+  (z  -  Zq)"1  Aj  +  A0(z)  (3.3-1) 
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where  AQ(z Q)  is  finite,  A^  /  0,  and  A^,  1  <  L  <  k 
are  constant  matrices.  If  Zq  =  oo  ,  (z  -  z^)"^  is  re¬ 
placed  by  z1,  1  <  i  <  k  .  All  of  AQ(z),  A^,  k^ 

are  uniquely  defined  by  their  construction  from  A(z)  . 

Definition  3.3-1.  If  A(z)  is  given  by  equation  3.3-l> 
then  k  is  the  order  of  the  pole  of  A(z)  at  z  =  zQ  . 

Definition  3.3-2.  A  complex  rational  matrix  is  said  to  bo 
reverse-hermltlan  if  A*(z)  =  A(l/  z)  (the  function  A  is 
symmetric  with  respect  to  the  unit  circle).  Hence,  on  the 
unit  circle,  z  -  elu),  A*(elui)  =  A(eiu))  and  A(eil")  is 
hermitian  in  the  ordinary  sense.  For  real  A(z)  ,  this 
condition  simplifies  to  A' (l/z)  =  A(z)  and  will  be  called 
reverse-symmetrical .  A  real  scalar  function  f(z)  satis¬ 
fying  f (l/z)  =  f (z)  is  also  called  reverse-symmetrical . 

It  is  most  convenient  for  typographical  reasons 

to  let 

A*(z )  s  A*(l/  z). 

This  notation  is  used  throughout  the  remainder  of  this 
paper.  Notice  that  A**(z)  =  A(z),  (A  B)#  =  B#  A*  . 

Definition  3.3-3.  A  rational  m  x  n  matrix  A(z)  is  said 
to  be  reverse-unitary  if 

A(z)  A*(z)  =  lm  . 

A  reverse-unitary  matrix  is  also  called  all-pass . 
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Definition  3. 3-4 .  A  matrix  A(z)  is  said  to  be  regular 
if  it  is  analytic  inside  the  unit  circle  |  z  j  <  1  .  A 
matrix  A(z)  is  said  to  be  Hurwltzlan  if  it  is  analytic 
inside  and  on  the  unit  circle  |  z  |  1  . 

3.32  Spectral  Decomposition 

The  decomposition  of  polynomial  matrices  that  is 
discussed  in  this  section  is  very  closely  related  to  that 
of  the  Spectral  Theorem  of  Linear  Algebra  (Hoffmai  and 
Kunze,  1961,  pp.  275-6)  which  is  stated  for  normal  opera¬ 
tors.  Thus  we  will  call  the  decomposition  theorem  the 
Spectral  Theorem. 

Before  stating  this  theorem  we  shall  investigate 
the  properties  of  the  latent  zeros  and  vectors  of  a  poly¬ 
nomial  matrix.  These  properties  will  account  for  the 
principle  restrictions  placed  upon  the  theorem. 

3.321  Latent  zeros  and  vectors 

Let  us  consider  the  n  x  n  square  polynomial 

matrix 

A  (z )  =  Ao  +  Alz  +  •••  +  Amz 

The  latent  zeros  of  A(z)  are  those  values  of 

z  =  zi  i  =  1,  ...,  p  (p  =  nm  if  A(z)  is  non-degenerate) 
for  which  Det  A(z)  =  0  .  Since  the  determinant  has  real 
coefficients,  complex  roots  may  only  occur  in  conjugate 
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pairs . 


Frazer,  et  al . ,  (pp.  61-65,  19^7)  prove  the 
following  properties  concerning  polynomial  matrices  at  the 
zero  positions  z^  : 

(a)  The  matrix  A(z^)  is  necessarily  singular. 
When  z^  * s  an  unrepeated  root,  A(z^)  has  rank 

r (  A(z^)  )  =  n  -  1  . 

(b)  When  A(z^)  has  rank  r(  A(z^)  )  =  n  -  q, 

at  least  q  of  the  roots  z^,  Zg,  z^  are  equal  to 

(c)  The  matrix  A(z^)  does  not  necessarily  have 
rank  n  -  q  when  z^  is  a  root  of  multiplicity  q  . 

(d)  When  A(z^)  has  rank  r(  A(z^)  )  =  n  -  1 

the  adjugate  Adj  A(z^)  has  unit  rank,  r(Adj  A(z^)  )  =  1 . 
Hence  it  is  expressible  as  a  product  of  the  form 

Adj  A(z^)  =  u^  v» 

where  u^  and  are  column  vectors  (called  latent 

vectors)  of  length  n  and  are  constants  appropriate  to  the 
selected  zero  z^  .  At  least  one  element  of  each  vector 
is  non-zero. 


Finally,  if  we  let 


Dm  A (z. )  = 


i~*n 


A  (z ) 


—  t-t  , 
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we  have 


(e)  When  A(z.)  has  rank  r(  A(z.)  )  =  n  -  q, 
where  q  >  1,  the  adjugate  matrix  Adj  A(z)  and  its 
derivatives  up  to  and  including  Dq~2  Adj  A(z^)  are  all 
null.  However,  the  matrix  Dq~  Adj  A(z^)  has  rank  q 
and  is  expressible  as  a  product  of  the  form 

Dq_1  Adj  A(zt)  =  jL,, 

where  and  are  n  x  q  matrices.  The  columns  of 

these  matrices  can  then  be  used  to  form  q  pairs  of  latent 
vectors  u^  and  . 

3.322  Spectral  Theorem 

It  is  frequently  convenient  to  introduce  the 
concept  of  2-term  operators  which  correspond  to  polynomial 
matrices  of  degree  1.  If  we  examine  one  of  these  2-term 
operators 

I  -  Uz  , 

we  see  that  it  is  closely  related  to  the  characteristic 
value  problem  that  is  usually  formulated  in  terms  of  A  : 

u  -  IX. 

Thus  we  may  apply  our  existing  knowledge  of  the  character¬ 
istic  zeros  and  vectors  of  constant  matrices  to  the  more 
general  case  of  polynomial  matrices.  This  approach  is  used 
in  the  spectral  theorem. 
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Theorem  3.3-1  (Spectral) .  Let  A(z)  be  an  n  x  n  real 
polynomial  matrix  of  rank  n  and  degree  m 

A(z)  =  AQ  +  A1z  +  ...  +  Amzm  . 

Then  A(z)  may  be  represented  as 

A (z)  =  C0(z)  (I  -  UlZ)  ...  (I  -  U^z) 

=  ^j(z)  •••  G^(z)  (3.3-2) 

or  as 

A  ( z )  =  (I  ®£z)  •••  “  4-^z)  ^q(z) 

=  \(z)  ...  ^(z)  C0(z)  (3.3-3) 

where  Gq(z)  and  ^q(z)  are  elementary.  If,  for  every 
zero  z^  of  multiplicity  q,  r(  A(z^)  )  =  n  -  q.^  where 
q/qx  <  l  . 

Proof.  (Claerbout  (personal  communication)  has  de /eloped 
a  similar  factorization.) 

First,  consider  equation  3.3-2.  Since 
j  A  B  j  =  |  A  |  |  B  |  ,  the  latent  roots  z i  i  =  1,  . . . ,  p 

of  A(z)  must  be  the  union  of  the  latent  roots  of 
G^,  G2,  ...,  G^  .  The  n  latent  vec^rs  of  A(z)  and 
G^(z)  are  given  by 

Adj  A  (z j)  =  Adj  (Gq  G1  ...  C^) 

=  Adj  Gt(Zj)  Adj  (Gq  G1  ...  C^1) 

“  ” j  “j  A^j  (G0  G^  ...  G^^ ^ )  j  =  1,  . <  < ,  n  . 
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(3.3-4) 


Thus  the  n  latent  vectors  3ij  are  the  same  for  A(z^) 
ai.u  .  Therefore,  if  we  determine  a  set  of  n 

zeros  of  A(z)  that  have  n  independent  latent  vectors, 
we  may  recombine  them  by  the  well-known  formula  (Frazer,  et 
al.,  pp.  66-68) 


"zn  0  ” 

1 

= 

(-l)  •  *(-u-n) 

0  " "  "z 

n 

(3.3-5) 

and 

"  1  ~  UtZ  * 


The  2-term  polynomial  0^(z)  is  a  factor  of  A(z).  For, 
if  we  substitute  the  matrix  into  the  polynomial  A(z) 

A(ut)  =  A0  +  Ax  ut  +  ...  +  Am  uj  (3.3-6) 


-  A0  +  A. 


Ut  Ut1 


A  u 
m 


-lyn  u-1 


l'*  l  1  t 


we  see  that  A(lt^)  =  0  identically.  Thus  G^(z)  =  I  -  U^z 
will  right-divide  A(z)  with  a  remainder  of  zero  (Frazer, 
et  al.,  19^7>  p.  60 )  and  therefore  is  a  factor  of  A(z). 

(Q.E.D.) 


The  factorization  is  continued  then  by  removing 
G^(z)  by  right  division,  determining  n  more  independent 
latent  vectors  and  constructing  a  second  2-term  wavelet 
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.  This  process  is  repeated  until  Gq(z),  an  elemen¬ 
tary  matrix,  remains. 


Alternately,  we  may  factor  on  the  basis  of  the 
latent  vectors.  Thus,  let  us  consider  equation  3 •  3“ 3 • 
Here  again  the  latent  zeros  of  A(z)  and  the  factors 
^(z)  are  the  same.  The  latent  vectors  of  A’  (z)  are 
given  by 

Adj  a1  (Zj)  =  Adj  ...  tr1  ?r0)' 

=  Adj  (I|(z  )  Adj  (cEt_1  ...  Cx  Z0)' 

-  Oj  uj  Adj  1?(.J  ...  fj  s0)'  (3.3-7) 

Thus  t\_.  ,,  latui  t  vectors  of  A'(z)  corrcr.poniin :  to 

z.,  i  =*  1,  2,  ...,  n  are  the  same  as  the  n  latent 
vectors  of  £  (z)  .  As  before,  if  we  choose  n  zeros 
such  that  the  associated  latent  vectors  are  independent 
then  they  may  be  recombined  as 


and  used  to  remove  Cl^(z)  from  A(z) 


1  _1 

(jLj )  •  •  •  (jLn ) 

(3.3-8) 
by  left  division. 


If  the  ze^s  and  vectors  are  independent  then 
the  choice  of  which  n  zeros  to  associate  with  each  2-term 
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Figure  3.3  -  1:  Spectral  decomposition  of  a  full  rank  polynomial  matrix  according 
to  Theorem  3.3-1* 
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-  o 
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-  O  -n 

'  *  1- 


z 

H 
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factor  is  arbitrary.  Altogether,  there  may  be  (nm)’./^n,^m 
different  factorizations.  Once  the  choice  is  made,  of 
course,  the  order  of  factorization  must  be  preserved  since 
multiplication  is  not  commutative,  in  general.  In  some 
instances  the  choice  of  zeros  must  be  made  under  certain 
restrictions  so  that  the  full  factorization  may  be  realized 
As  indicated  above,  this  restriction  consists  of  choosing 
the  zeros  so  that  the  latent  vectors  or  are  inde¬ 

pendent  for  each  set  of  n  vectors.  Such  a  choice  may 
always  be  made  if  for  every  zero  of  multiplicity 

q^,  r(  A(z^)  )  =  n  -  q^  where  q/q^  £  t  .  Q.E.D. 

The  details  of  this  factorization  are  illustrated 
in  Figure  3«3“1.  The  right  half  shows  the  decomposition 
in  terms  of  and  the  left  half  shows  the  decomposition 

in  terms  of  . 

In  general  It  is  not  necessary  to  go  through  the 
intermediate  steps  of  forming  2-term  factors  to  construct 
a  polynomial  matrix  from  its  latent  roots  and  vectors. 

This  more  direct  approach  is  the  subject  of  the  Spectral 
Corollary. 

Corollary  3.3-1  (Spectral) .  Let  A(z)  be  an  n  x  n  real 
polynomial  matrix  of  rank  n  and  degree  m 

A(z)  =  kQ  +  Ax  z  +  . . .  +  Am  z 


41 


■  1  *  i  p^guggpp 


Then  A(z)  is  completely  described  by 

a)  the  elementary  matrix  riQ(z) 

b)  the  latent  zeros  z^  i  =  1,  2,  ...,  p,  and 

c)  the  p  (p  <  m  n)  corresponding  latent 

vectors  or  v.  if  for  every  root  z^ 

of  multiplicity  q,  r(  A(z^)  )  =  n  -  q  . 

Notice  that  this  corollary  is  not  so  general  as 
the  spectral  theorem  in  its  treatment  of  multiple  zeros 
with  identical  latent  vectors. 

Proof  (Suggested  by  Quenoullle,  1957*  PP«  5-25)  Let  us 
firzi  consider  the  case  for  which  A(z)  is  non-degenerate, 
l.e.  that  GQ(z)  is  a  constant  non-singular  matrix  and  the 
degree  of  the  determinant,  Det  A(z),  is  p  =  m  n  . 

Consider  the  factored  form  of  A(z) 

A(z)  -  Aq  (I  -  Ux  z  -  ...  -  Um  zra) 

-  A0  C(z)  (3.3-9) 

where  .  Then,  if  we  inquire  about  the  solu- 

tions  to  the  equation 

u  -  U1  u  z  -  ...  -  Um  u  zm  =  0,  (3.3-10) 

we  see  that  solutions  are  possible  only  if  the  determinant 

I  -  U:  z  -  ...  -  Um  zm  =  0  (3.3-11) 

is  zero.  It  is  zero  at  the  p  locations  z^  i  =  3,  ...,  p 
which  are  the  latent  zeros  of  G(z)  and  consequently  of 


U2 


A(z)  .  Therefore  the  solution  vectors  of  equation  3.3-10 
are  the  latent  vectors  of  r(z).  Since 

Adj  A  ( z j  )  -  Adj  (AQ  G(Z<t)) 

=  Adj  r( zt)  Adj  Aq 

=  u^  v|  Adj  Aq  ,  (3.3-12) 

the  solution  vectors  are  also  latent  vectors  of  A(z)  . 

Nov/  join  the  latent  vectors  and  latent  zeros  into  the 
modal  matrix 


Clearly  we  can  solve  this  set  of  simultaneous  equations  for 
U2,  . . .  if  the  columns  of  li  Z  are  independent. 
This  does  not  occur  when  a  zero  of  multiplicity  q  has 
fewer  than  q  independent  latent  vectors,  that  is,  if  a 
multiple  root  has  Identical  latent  vectors. 


Alternately,  we  may  choose  to  use  the  vectors 
_v^  for  the  reconstruction.  For  this  case  we  would  factor 
A(z)  as 


A3 


A(z) 


=  (I  -  V1  z  ...  um  z  a)  A0 

=  (z)  Aq  (3.3-lb) 

Where  V\  =  -  A^  .  Since  the  latent  zeros  of  £(z) 

are  the  same  as  for  A(z),  and  since  the  latent  vectors 
of  A'  are 

Adj  A 1 (z^)  =  Adj  (?(z,)  Aq) ' 

=  Adj  ?•»  (z,)  Adj  A^ 

=  v  u'  Adj  A^  (3.3-16) 

we  may  reconstruct  the  matrix  G(z)  by  the  simultaneous 
equations 

»|VZr  ,.,tv;i.z"  =  1<  (3.3-17) 

where 

lr  =  (v2)  (v2)  ...  (vp)  .  (3.3-18) 

mm  “ 

The  same  restraints  hold,  here  as  held  for  the  u  vectors. 

If  Gq(z),  the  elementary  matrix  multiplier,  is 
not  constant,  then  the  number  of  zeros  is  p  <  mn  .  We  can, 
however,  proceed  as  above  to  find  G(z)  and  then  deter¬ 
mine  Gq(z)  by  the  formula 

Oq(z)  =  A (z )  G"1(z) 

for  the  factorization  ir.  terms  of  the  u  vectors,  or  by 
the  f ormu  j.a 

Cl0(z)  =  ft"1  (z )  A  (z ) 

<4b 


for  the  factorization  in  terms  of  the  _v  vectors. 

Q.E.D. 

These  factorizations  are  illustrated  in  Figure 
3.3“2.  The  boxes  enclosed  in  dottea  lines  represent  com¬ 
pletely  equivalent  representations  of  the  matrix. 


Example  3.3-1 .  (after  Claerbout,  personal  communication) 


A (zC )  - 


Consider  the  polynomial  matrix 


2  -  20z  +  50z2  -1  +  9z  -  20z2 
l4z  -  58z2  1  -  llz  +  28z2 


The  determinant  is 

A(z)  =  2  -  28z  +  l42z2  -  308z3  +  240zi* 

=  2  (1  -  2z)  (1  -  3z)  (1  -  Hz)  (1  -  5z) 

The  adjugate  matrix  is 


Adj  a (  )  = 


1  -  llz  +  28z2  1  -  9z  +  20z2 

-  I4z  +  58z2  2  -  20z  +  50z2 


(3.3-19) 


(3.3-20) 

Substituting  z^  =  0.5  into  Adj  A(z)  gives  Lhe  latent 


vectors  u^  and  v^ 


Adj  A (0.5)  = 


2.5  1.5 

7.5  4.5 


0.5 


1 

3 


[5  2] 


Likewise,  if  we  substitute  for  the  other  roots  we  will  find 


all  the  latent  vectors: 


Zero 

Latent 

vector 

1/2 

(1, 

3)' 

1/3 

(1* 

4) ' 

1/4 

(o, 

1)' 

1/5 

(1. 

6) ' 

u  Latent  vector  v 

(5,  3)' 

(2,  1)’ 

(1.  1)' 

(1.  0)' 


Now  let  us  follow  the  reconstructions  considered  in  the 
Spectral  Corollary. 


First,  we  may  pre-dlvlde  by  A„  to  find 


A(z) 


2  -l" 

/ 

"3l“ 

‘  4  -4* 

( 1  - 

z  - 

0  1 

J 

\ 

-14  11 

_  — 

_58  -28 

2 

2 


(3.3-21) 


Set  up  the  matrices  U  and  7  and  substitute  Into  the 


the  transpose  of  equation  3.3-14 


1/2 

3/2 

iA 

3  A 

1/3 

4/3 

1/9 

4/9 

0 

1 A 

0 

1/16 

1/5 

6/5 

1/25 

6/25 

1  3 

1  4 

0  1 
1  6 

(3.3-22) 


and  obtain  the  values  of  and  given  above  In 

equation  3.3-21. 


Similarly,  we  may  post-divide  A(z)  by  Aq  to  find 
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A  (z)  = 


I  - 


10  1 
L-7  4. 


-2S 


■5 


2  9  1 


-1 


0 


I  -  V1z  -  V2 z‘ 


‘0 


(3.3-23) 


We 

set  up 

the 

matri 

ces  If 

equation 

3.3- 

•17 

"5/2 

3/? 

5/4 

3/4  “ 

2/3 

1/3 

2/9 

1/9 

1/9 

1/4 

1/16 

1/16 

1/5 

0 

1/2^ 

0 

to  obtain  the  values  of 
equation  3  3-23* 


and  7  and  substitute  Into 


~  L, 

> 

3' 

vi 

2 

1 

1 

1 

V2 

1 

0 

Mb  m 

M. 

_ 

(3.3-24) 

and  V2  given  above  in 


Equations  3.3-21  or  3.3~23  may  also  be 
reconstructed  using  the  algorithm  of  the  Spectral  Theorem. 
We  will  illustrate  the  process  for  only  the  u  latent 
vectors.  Recall  that  this  factorization  is  in  the  form 

A(z)  =  AQ  (I  -  U1  z)  (I  -  U2  z) 

=  Aq  G1(z)  Gp (z )  . 

We  arbitrarily  choose  the  zeros  z  =  1/4,  1/5  to  obtain 


"2  - 

1 

o 

h-1 

_ i 

i 

o 

1 _ 

i - 

r*i 

o 

i _ 

.1  6- 

-0  5. 

— i 

o 

1 — 1 

_ 1 

5  0 
6  4 


*4  r ;• 


and 


I  - 


5  0 
6  4 


If  we  now  post  divide  G( z)  (from  equation  3.3-21)  we  find 


c. 


I  - 


-2  1 

-20  7 


which  has  latent  vectors 


zero 


1/2 

1/3 


latent  vector  u 

(1,  -4) ' 

(1.  -5)' 


Of  course,  this  is  not  the  only  possible  factori¬ 
zation.  Altogether,  there  are 

Ml.  =  24  _  6 

(ni)m  4 

different  representations.  Using  the  method  illustrated 
above  we  find  that 
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Example  3.3-2 .  (Multiple  roots) 


Consider  the  polynomial  matrix 


A(z)  = 


3  -  6  z  +  3  z‘ 
1  -  2  z  +  z^ 


2  -i 


1  -  4  z  +  z 

2  -  8  z  +  7  z‘ 


Using  the  standard  factorization  techniques  we  find  that 
this  has  latent  zeros  and  vectors 
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zero 


latent  vector  u 


1/2 

(i.  i) 

1/2 

(i.  i) 

1 

(1.  0) 

1 

(1,  0) 

If  we  set  up  the  reconstruction  equation  3.3-1^ 


1/2 

1/2 

i  A 

iA 

1/2 

1/2 

i  A 

i  A 

1 

0 

l 

0 

1 

0 

i 

0 

1  1 
1  1 
1  0 
1  0 


we  find  that  the  left  hand  side  is  singular.  However,  if 
we  use  the  algorithm  outlined  for  the  Spectral  Theorem  and 
use  the  zeros  1/2  and  1  for  each  of  the  2-term  factors, 
we  find 


A(z)  = 


Thus,  this  approach  is  slightly  more  general. 


3.33  Invertibility 

Definition  3.3-5.  A  one-sided  matrix-valued  wavelet  A(z) 
is  said  to  be  Invertible  if  there  exists  a  one-sided  left- 
or  right-inverse  wavelet  A  ^(z)  . 

Let  us  consider  only  square  matrix-valued  wavelets. 
The  inverse  of  such  a  wavelet  is  given  by 

A-1 (z )  =  AdJ  A(z) 

Det  A  (z)  . 


The  condition  for  invertibility  is  that  the  determinant 
of  A(z)  has  a  stable  inverse.  This  condition  when  applied 
to  the  determinant  of  a  finite  wavelet  is  exactly  the  same 
as  that  applied  to  the  scalar  wavelet.  That  is,  the  zeros 
of  the  determinant  of  A(z)  must  be  outside  the  unit 
circle  in  the  z  plane  (see  section  3.22). 

3.3^  Smith-McMillan  Canonical  Form 

This  canonical  form  for  rational  matrices  in¬ 
volves  the  terms  contained  in  the  determinant  and  the  rank 
of  the  matrix.  It  is  the  subject  of  the  classical  Smith- 
McMillan  Theorem  (Gantmacher,  p.  134,  1959  and  McMillan, 
p.  581,  1952). 


Theorem  3.3-2  (Smith-McMillan) .  Let  A(z)  be  an  m  x  n 
complex  rational  matrix  of  normal  rank  r  .  Then  there 
exist  two  elementary  polynomial  matrices  C(z)  and  F(z) 
of  orders  m  x  r  and  r  x  n,  respectively,  such  that 


A(z) 


C (z )  diae 


^(z) 

*2  (z) 


CDF 


m2(z) 

’^(z) 


”r(z)  " 
'''r  (z)  _ 


F(z) 


(3.3-25) 


where 


a)  ^(z)  and  '!,^(z)  are  relatlvely  prime 
polynomials  with  unit  leading  coefficients,  1  <  k  <  r; 


5  2 


j  in  1.  .1 


b)  Each  ro^(z)  divides  rDj<+i(z)> 

1  <  k  <  r  -  1,  and  each  ^(z)  divides 

2  <  i  <  r; 

c)  The  diagonal  matrix  D(z)  appearing  in 
equation  3.3-25  is,  subject  to  a)  and  b),  uniquely 
determined  by  A(z)  .  It  is,  in  fact  canonic; 

d)  If  A (z)  is  real,  the  n’s,  *'s,  C(z), 
and  F(z)  may  also  be  chosen  real; 

e)  The  finite  point  z  =  Zq  is  a  pole  of  A(z) 
of  order  k  if  and  only  if  it  is  a  zero  of  '^(z)  of 
order  k  . 

f)  The  order  of  z  =  co  as  a  pole  of  A(z)  is 
the  same  as  the  order  of  1/z  =  0  as  a  pole  of  A(l/z)  . 

A  rational  matrix  is  said  to  be  Smlth-McMillan 
canonic  if  it  is  square,  non-singular  and  diagonal  with 
properties  a)  and  b)  listed  above.  The  rational 
functions  ^r/*r  are  generalized 

invariant  factors  of  A(z)  .  Clearly,  since  C  and  D 
are  elementary,  Det  A(z)  =  Det  D(z)  .  A  set  of  polynomials 
are  said  to  be  relatively  prime  if  their  largest  conmon 
denominator  is  1. 

Frazer,  et  al .,  (pp.  87~92,  19^7)  or  Gantmacher 
(pp.  13^-139j  1959)  show  in  detail  the  technique  for  the 
reduction  of  a  matrix  to  canonical  form.  The  method  em¬ 
ployed  is  reminiscent  of  the  elimination  methods  for  in- 
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vertinc  a  matrix.  We  will  Illustrate  it  with  a  polynomial 
matrix  example. 

Example  3.3-3.  Reduction  of  a  polynomial  matrix  to  can¬ 
onical  form. 


Let  us  consider  the  matrix 


A  (z) 


2  +  z 
1 


z 

6  +  z 


a)  Put  a  one  in  the  first  diagonal  position.  This  is 
accomplished  by  multiplying  A  on  the  left  by 


*  ■  [; :] 


to  obtain 


1  6  +  zl 

_2  +  z  z  J  . 


b)  Reduce  the  other  terras  in  the  first  row  and  column  to 
zero.  This  accomplished  by  multiplying  on  the  left  by 


to  obtain 


z) 


0 

1 


6  +  z 

-z2  -  7z  -  12 


and  multiplying  on  the  right  by 


BLANK  PAGE 


0^  “  1  -  (6  +  z  )"J 

Lo  -i  J 

to  obtain 

1  °  1 

0  z2  +  7z  +  12J 

c)  Now,  if  we  let 

;] 

F  (z)  .  O'1  -  r  1  -6  -  z 

Lo  -i 

we  obtain  the  Smith-McMlllan  canonical  form 
A  (z )  *  CDF 

where 

1  °  i 

0  (z  +  3)(z  +  4)J 

Clearly  C  and  F  are  elementary  matrices  and 
D  is  canonic.  That  is,  1  divides  (z  +  3)(z  +  4)  and 

Det  D(z)  =  Det  A(z)  . 
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3.35  Robinson  Canonical  Form  and  All-Pass  Systems 

Theorem  3.3*3  (Robinson  Canonical  Form) .  Any  full-rank 
wavelet  A(z)  can  always  be  uniquely  represented  by  the 
Robinson  canonical  form 

A(z)  -  AQ(z)  P(z) 

where  A^(z)  is  invertible  and  P(z)  is  regular  reverse¬ 
unitary  (l.e.  all-pass).  More  generally,  if  A(z)  is  an 
n  x  m  matrix  and  has  rank  r  <  n,  m,  then  Its  canonical 
form  becomes 


where  P(z)  is  a  regular  reverse-unitary  m  x  m  matrix. 

Matrix  re verse -unitary  (all-pass)  systems  have 
similar  properties  to  scalar  all  pass  systems.  V,’o  shall 
state  several  theorems  concerning  them  now  (Robinson. 

1962,  and  You la,  1962). 

Theorem  3.3-^.  An  n  x  m  matrix  P(z)  of  rank  r  is 
analytic  in  the  entire  z  plane  together  with  its  Inverse 
(either  left,  right,  or  both)  if  and  only  if  it  Is  an 
elementary  polynomial  matrix. 

Proof:  The  "if”  part  is  obvious.  According  to  the 
Smith-McMillan  Theorem  (3.3“2),  the  analyticity  of  P(z) 
for  all  z  implies  that  all  of  the  denominator  terms, 

of  the  canonical  form  be  constant.  Now  the  existence 


of  a  left  or  right  inverse  implies  that  either  n  =  r  or 
m  *  r,  respectively.  The  canonic  form  for  p”*(z)  is 


dlag 


*r(z)  ♦r(z) 


*1  (z) 


l«rM  Vl^^  ^i(z)  J 


The  analytlclty  of  P“1(z)  in  the  entire  plane  implies 
that  cp^,  t  ■  1,  . ..,  r  is  constant.  Therefore  P(z) 
is  the  product  of  three  elementary  polynomial  matrices, 
of  rank  r  .  Q.E.D. 


Theorem  3.3-5.  A  reverse -unitary  rational  matrix  is 
bounded  on  the  unit  circle. 


Proof:  Suppose  P(z)  is  m  x  n  and  P(z)  P#(z)  -  1  . 

Thus  P(eltw)  P  (e1<w)  -  ln,  and,  writing  out  the  diagonal 
elements  ln  expanded  form, 

^  |  (?)rk^el<W)  |  2  *  1  (^  *  1*  2,  ...,  n)  . 

|(p)rk(ela,)|  <1  (r  -  1,  2,  ...,  mj 

k*  1,  2,  ...,  n), 

for  all  iH  .  Q.E.D. 

Theorem  3.3-6.  The  only  regular  reverse-unitary  matrices 
P(z)  with  regular  inverses  are  constant  v  .tary  matrices 
(trivial  all-pass  systems).  If  P(z)  is  real  it  is  real- 
orthogonal  . 

Proof:  Suppose  P(z)  P#(z)  -  1R,  say,  where  P(z)  is  a 
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regular  n  x  m  reverse-unitary  matrix.  The  analyticity 
of  its  right  inverse  inside  the  unit  circle  implies  that 
of  "P(l/  z)  in  the  same  region  and  therefore  that  of 
P(z)  outside  the  unit  circle  including  infinity.  Now  the 
poles  of  P(z)  are  the  complex  conjugates  of  those  of 
P(z)  .  Hence  P(z)  is  analytic  in  the  entire  z  plane 
and  bounded  at  infinity.  By  Liouville's  Theorem  it  must 
be  a  constant  unitary  matrix.  If  P(z)  is  real  it  must 
be  real  orthogonal  by  definition.  Q.E.D. 


3.36  Delay 

The  delay  of  a  one-sided  matrix-valued  wavelet 
is  a  measure  of  how  the  operator  redistributes  the 
energy  of  an  input  process.  It  is  defined  in  terms  of 
the  partial  energy 


i 


The  following  theorem  is  a  discrete  analog  of  a 
theorem  given  by  Robinson  (pp. 83-88,  1962).  Since  his 
proof  is  rather  long  and  involved,  it  will  not-  be  repeated 
here. 


Theorem  3*3-7  (Minimum  Delay) .  The  delays  of  the  set  of 
wavelets  Ai(z)  which  have  the  same  invertible  Robinson 
canonical  form  Aq(z)  are  greater  than  or  equal  to  the 
delay  of  AQ(z),  Equality  holds  if  and  only  if  the  all- 
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3>3? 


pass  system  P(z)  is  trivial.  That  is,  the  partial  ener¬ 
gies  obey  the  relation 

j0  tr  i  j0  <AO>J<A0>j‘  *11  * 

where  J  is  the  time  index. 

3.37  Phase 

As  in  the  scalar  case,  the  Fourier  transform  of 
the  operator  A(z)  is  determined  by  restricting  our  atten¬ 
tion  to  z  =  e  .  We  may  proceed  to  express  each  poly¬ 
nomial  element  of  the  matrix  in  terms  of  an  amplitude 
characteristic  and  a  phase  characteristic.  The  question 
then  arises  whether  there  are  any  simple  measures  of  this 
phase  matrix,  other  than  the  determinant  of  the  polynomial 
matrix,  which  would  correspond  to  invertibility.  That  is, 
can  we  formulate  a  minimum  phase  theorem  for  matrix-valued 
wavelets  ? 

An  empirical  investigation  was  made  of  this  ques¬ 
tion  which  gave  negative  answers  for  all  measures  tried. 
These  measures  included  1)  the  trace  of  the  phase-lag 
matrix,  2)  the  phase-lag  of  the  trace  of  the  polynomial 
matrix,  and  3)  the  norm  of  the  phase-lag  matrix  as  a 
function  of  .  In  every  case  counter-examples  could  be 
found  for  which  the  wavelet  was  invertible  but  the  measure 
tried  did  not  give  a  minimum. 
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3.4 


Another  approach  which  might  prove  more  success¬ 
ful  would  be  to  define  a  matrix  amplitude  characteristic 
Q(ui)  and  a  matrix  phase  characteristic  such  that 

*  /  ~i«>\  n  t  \  -i£(u>) 

A  (e  )  =  G(uj)  e  '  ' 

This  has  not  beer.  Investigated.  However,  should  it  prove 
to  have  a  minimum  phase  property  associated  with  inverti- 
bility,  this  measure  would  have  limited  application  because 
of  the  difficulty  of  computation  and  cognition  of  such 
characteristics . 

3*^  Multi -Dimensional  Wavelets 

Our  treatment  of  multi-dimensional  wavelets  will 
be  brief  on  two  accounts.  First,  there  is  no  general 
factorization  available  for  multi-dimensional  polynomials; 
and  second,  in  almost  all  problems  with  which  we  are  con¬ 
cerned  the  multi-dimensional  process  can  be  mapped  into  an 
equivalent  higher  ordered  one -dimensional  matrix-valued 
process. 


The  absence  of  any  factorization  can  be  illus¬ 
trated  by  attempting  to  factor  the  polynomial 


aoo 

+ 

aoi 

X 

+ 

a02 

xc 

+  aio 

y 

+ 

all 

xy 

+ 

al2 

x2y 

2 

2 

2  2 

+  a20 

y 

+ 

a21 

xy 

+ 

a22 

x  y 

60 


■ys* 


If  this  were  factorable,  we  should  be  able  to  find  two 
polynomials. 


b(x,y)  = 

b00  +  b01x 

and 

'"(x,y)  =  cQ0  +  cQ1x 

b10y  +  bllxy 

c10y  +  cUxy 

such  that 

b  c  =  a  .  But 

a  has 

9  degrees  of  freedom  and 

b  and  c 

combined  have 

only  8. 

Thus,  unless  there  are 

special  relationships  between  the  elements  a. a(x,y) 
is  cnfactorable . 

3.41  Invertibility 

Definition  3.4-1.  A  scalar  multi-dimensional  wavelet 
a(z)  is  said  to  be  invertible  about  its  origin  z  =  0 
if  there  exists  a  wavelet  a”1^)  such  that 
a(z)  a"1^)  »  1  . 

The  condition  for  making  an  expansion  about 
z  =  0  of  a  2(z)  is  that  a(z)  does  not  go  to  zero 
inside  the  unit  hypercircle 

|  Z  Zj  ...  znj  .  1  . 


3.42  Phase 

Perhaps  the  simplest  measure  of  invertibility 
involves  the  phase-lag  of  the  wavelet.  The  multi¬ 
dimensional  Fourier  transform  is  found  by  restricting  z 
to  the  unit  hypercircle.  Thus 
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-itD  -iui^  -ian 

a  (zjzjj  zn)  ^  a(e  >  e  J  e  ) 

and  by  finding  the  polar  representation  of  this  we  can  find 
a  multi-dimensional  phase-lag  characteristic, 

-m(uj,  oj^,  ...f  ui^)  . 

Theorem  2.4- a  (Minimum  Phase) .  If  the  phase-lag  character¬ 
istic  -o(w)  for  the  wavelet  a(|J  is  the  same  for  all 
a)  «  n  or  0  i  =  1,  . ..,  n,  then  the  wavelet  a(z)  is 
invertible . 


Figvjre  3.',-l  illustrates  two  phase-lag  plots  for 
two-dimensional  wavelets.  The  variable  z^  corresponds 
to  the  phase-variable  .  Notice  that  when  a  hyper¬ 
surface  a (2.)  »  0  cuts  across  the  unit  hypercircle,  the 
phase  is  discontinuous  along  the  intersection.  This  is 
analogous  io  the  case  of  a  zero  on  the  unit  circle  for 
one-dimensional  scalar  wavelets. 
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Invertible  wavelet 


Figure  3.4  -  Is  two  dimensional  phene- lag  characteristics 
for  two  wavelets.  The  phase* lag  is  discontinuous 
across  the  zero  hyper-surface  in  the  second  case. 
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3.43  Mapping  into  one-dimensional  representation 

Much  of  the  algebra  of  multi-dimensional  operators 
and  autocorrelations  represents  a  special  case  of  the  gen¬ 
eral  matrix- valued,  one-dimensicnal  algebra.  For  this 
reason,  we  seek  to  map  multi-dimens tonal  convolution  into 
a  matrix-valued  notation  rather  than  to  develop  the  algebra 
In  multi -dimensional  notation.  Thus,  this  section  will 
give  an  extensive  account  of  a  mapping  from  multi-dimensional 
notations  to  one-dimensional  notation. 

As  pointed  out  in  Chapter  2,  this  mapping 
necessarily  assumes  a  preferred,  or  time-like,  direction. 

It  is  this  dimension  that  remains  undisturbed  after  the 
mapping.  Thus,  rather  than  thinking  of  a  multi¬ 
dimensional  wavelet  as  a  lump  in  multi-dimensional  space, 
we  may  visualize  it  as  a  set  of  time-wavelets  associated 
with  various  spatial  positions.  Then  we  take  the  logical 
step  of  placing  these  time-wavelets  into  a  vector  repre¬ 
sentation.  This  process  is  illustrated  for  a  three- 
dimensional  wavelet  in  Figure  3.4-2.  Notice  that  before 
the  vector  representation  can  be  accomplished,  we  must 
make  some  arbitrary  ordering  of  the  spatial  points. 

Now  let  us  consider  convolution.  If  the  opera¬ 
tor  wavelet,  a,  is  mapped  into  a  vector  of  time  wavelets, 
a,  and  the  output,  y,  is  mapped  into  a  similar  vector, 

£»  then  the  input,  x,  must  be  mapped  into  a  matrix,  X  . 
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a^z)  a4(z)  a5(z) 


Spread  these  time-wavelets  into  a 
vector  representation 


a(z)  *  ^a^(z)*  Sg( z)*  a^(z)»  a4(z)«  t^(z) j 


FInjre  3.M  -  2:  Mapping  of  a  multi-dimensional  wavelet 
into  vector  notation. 
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This  process  is  illustrated  in  Figure  3.4-3.  Each  column 
of  X  represents  the  configuration  for  the  dot  product 
for  one  spatial  lag  of  the  convolution.  We  can  think  of 
this  mapping  for  each  column  as  the  superposition  of  the 
spatially  reversed  a  grid  onto  the  x  grid  at  some  lag. 
The  lag  for  a  particular  column  corresponds  to  the  order¬ 
ing  of  the  output  grid. 

Let  us  now  put  the  methods  discussed  above  on 
a  formal  basis.  Since  there  are  no  well-defined  operators 
for  this  mapping,  we  will  use  short  mnemonic  words  to 
represent  each  operator.  These  operators  will  be  used 
only  to  define  such  mappings  as  described  above. 

In  nearly  all  of  these  discussions  the  scalar 
elements  may  be  replaced  by  matrices,  however,  to  avoid 
undue  confusion  we  will  make  the  definitions  in  terms  of 
scalar  quantities. 

REV  -  Reversing  operator 

REV  reverses  the  positive  sense  of  all  dimensions 
of  a  process: 

REV(X1)  «  X_, 

REV(X(z) )  =  X(l/z) 

where  we  define  -_i  =  (-i,  -i^,  •••>  -1^) 

\/7^  ■  (l/Z,  l/Z^i  •••#  l/Zjj)  • 
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3.43 


at,i1,i2  xt,i1#i2 


t,iri2 


•  »  convolution 


_ L 


z-fransform  in  Preferred  Direction 
and  order  the  indices. 

- T - 


ai(z)  a3(z)  #  x2(z)  x3(z) 

a2(z)  a4(z)  a5(z)  Xg(z) 


y^*)  y4(z)  y?(z) 

y2(*)  y5(*)  yqM  y10(2) 

y3(z)  y6(z)  y9(z) 


F 


Nap  wavelets  into  matrix  notation 
(a^  represents  ai(z)»  etc.) 

T 


2#  *2'  *51 


[ 


x1  Xg  0  x3  0  0  0  0  0  0 

0xlx20  x3  00000 

OOOx^XgOXjOOO 
0  0  0  0  Xj  Xg  0  x3  0  0 

0  0  0  0  0  0  0  Xl  Xg  x3 

yj  y2  y3  y4  y5  y6  y7  y8  y9  y10J 


Figure  3.4  -  3r  Mcppin^  of  multi-dimensional  convolution 
into  matrix  representation. 
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SHIFT  -  Origin  shifting  operator 

SHIFT  alters  the  origin  of  a  process  by  adding 
a  value  to  each  index  of  the  process: 

SHIFT  (J,  X^)  -  *1  +  j 

SHIFT  Q,  X(z))  .  X(z)  zJ  z^1  ...  zJN 

WINDOW  -  Window  operator 

WINDOW  isolates  a  portion  of  a  process,  Y, 
by  superimposing  the  grid  of  a  process  X  onto  the  grid 
of  Y  .  The  indexing  of  the  new  process  is  that  of  the 
window  grid  X  •  We  assume  that  Y  has  zeros  wherever  X 
extends  beyond  the  defined  limits  of  Y  . 


WINDOW 

(X.  X) 

*  z 

For  example  consider 

the  2-dimensional  process 

xi  - 

x0,-l 

x0,0 

x0,l 

*1.-1 

Xl,0 

Xl.l 

n 

>r' 

y-l,0 

y-l,l 

y0,0 

yo,i 

yi,o 

yl,l 

Then  the  WINDOW  operator  behaves  as 
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Z1  - 

WINDOW 

(X,  Y) 

- 

0 

y0,0 

y0,l 

0 

*1,0 

yl.l 

zo,-i 

Z0,Q 

Z0,1 

*1.-1 

*1.0 

*1.1 

ORDER  -  Ordering  operator 

ORDER 

converts 

a  vector  of  1  »  (t,  1^,  ...,  i^) 

into  another  vector  (t,J)  such 

that  J  takes  on  a 

unique  value  for 

each  of 

the  grid  positions  (i^,  .,,,  i^) 

of  a  finite  process: 

ORDER 

<xt 

*  • 

t ) 

-  xt.j 

The  actual  process  used  to  select  the  order  of 
enumeration  is  entirely  arbitrary  and  need  not  be  speci¬ 
fied  until  a  specific  application  is  made. 


MAPI  -  Mapping  operator 

MAPI  maps  a  multi-dimensional  process  into  a 
vector-valued  process.  Consider  a  multi-dimensional 
process  A*.  .  ■  ORDER  (A*  * )  1  <  J  <  N  then 

mpi  *  At 

■[At,l*  *t,2'  At.N  ]  * 

MAP2  -  Mapping  operator 

MAP2  corresponds  to  the  matrix  mapping  of  X(z) 
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that  was  made  in  Figure  It  is  defined  in  terms  of 

the  operators  above. 


X  = 

MAP2  (a,x) 

s 

MAPI  ^  ORDER  (Xt^k) j 

A 

where  Xt  k  = 

MAPI  ^ ORDER  ^WINDOW  ^  SHIFT (k,  a),  REV(x|^  ' 

A 

That  is,  X^  k 

represents  the  columns  of  the  matrix  X 

located  in  an  array  like  that  of  the  convolution  of  a 

A 

and  x  .  The  individual  columns  of  k  are  formed  by 


shifting  the  grid  of  a  by  an  amount  k,  superimposing  it 
on  the  spatial  reverse  of  x,  and  then  ordering  and  map¬ 
ping  this  intersection  according  to  the  Indices  of  a  . 

Now,  in  terms  of  the  z- transforms 

a(z)  x(z)  =  y(z) 

corresponds  exactly  to 

MAPI  (a(z) )  MAP2  (a,  x(z))  =  MAPI  (y (z) ) 

where  a(z)  represents  the  z-transform  of  a  in  the 
preferred  direction  only. 
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4.  FACTORIZATION  OF  AUTOCORRELATIONS 

The  operation  of  autocorrelation  is  generally 
defined  as  the  expected  value  of  the  cross-product  of  a 
process  with  Itself  as  a  function  of  time  and  spatial  lags. 
It  has  the  very  useful  property  of  removing  all  phase  in¬ 
formation  from  a  stochastic  process.  If  a  time-series  may 
be  characterized  as  the  convolution  of  a  white  light  process 
with  a  wavelet,  then  the  autocorrelation  of  the  process 
isolates  the  amplitude  properties  of  the  wavelet.  This  is 
because  the  autocorrelation  of  white  light  is  zero  except 
for  a  pulse  at  lag  zero.  These  properties  of  stochastic 
processes  have  been  treated  by  many  authors  and  from  many 
different  viewpoints.  Some  of  the  salient  works  include 
Riesz  (1907  and  1952),  FeJer  (1916),  Kolmogorov  (1941a  and 
b),  Karhuenen  (1947  and  1$49),  and  Szego  (1959).  Wold 
(1938)  stated  the  decompositlonal  properties  In  terms  of 
stochastic  time  series  as  follows: 

Theorem  4.1-1  (Wold  Decomposition) .  Any  stationary  process 
Xt  can  be  uniquely  represented  as  the  sum  of  two  mutually 
uncorrelated  process  *  Ut  +  Vt,  where  Ut  Is  deter¬ 
ministic,  and  Vt  is  the  convolution  of  a  one-sided  wave¬ 
let  with  a  stationary  white-light  process. 

Robinson  (1962)  and  Wiener  and  Masanl  (1957)  have 
extended  this  theorem  to  specify  a  particular  decomposition 
in  terms  of  an  Invertible  wavelet. 
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4.1 

With  this  brief  discussion  of  the  motivation 
of  autocorrelations  (for  more  detailed  discussions  see 
Wiener,  1949;  Whittle.  1954;  Wiener  and  Masani,  19b7  and 
19*33;  and  A, istev .  1964)  we  will  go  directly  to  a  discussion 
of  their  properties  and  factorizations.  In  general,  most 
of  the  factorizations  are  made  in  terms  of  correlations 
of  finite  length;  however,  some  of  the  cases  are  easily 
extendible  to  infinite  lengths. 

4.1  One-Dimensional  Scalar  Autocorrelations 

The  theory  o'"  one-dimensional  scalar  autocorrela¬ 
tions  is  well  known.  Thus  we  need  only  state  results  in 
this  section  for  the  purpose  of  giving  an  intuitive  intro¬ 
duction  to  the  following  sections. 

Let  r(z)  represent  a  real  autocorrelation  of 
length  m+  m  +  1 

r(z)  -  r-m  z‘m  +  ...  +  r-i  z'1  +  r0  +  rx  z  +  ...  +  rm  z™ 
then 

a)  r(z)  is  reverse-symmetric,  that  is 

r(z)  =  r#(z) 

-  r(l/z )  . 

b)  r(e’’i',))  is  non-negative,  that  is,  the 
cosine  transform  of  the  autocorrelation  is 
non-negative . 
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c)  ro  —  ri  with  equality  holding  only  if  the 
input  process  Is  periodic,  i.e.  deterministic. 

d)  The  real  frecuency  zeros,  that  is,  the  zeros 

on  the  unit  circle,  j  z  |  *  1,  are  of  even 

multiplicity. 

e)  For  every  zero  zi  of  r(z)  inside  the  unit 
circle,  there  is  a  corresponding  zero  1/z^ 
outside  the  unit  circle. 

It  is  interesting  to  note  that  since  the  cosine 
transforms  of  autocorrelation  functions  and  of  minimum  de¬ 
lay  wavelets  are  both  non-negative  (see  Section  3-25) >  the 
center  point  and  right  half  of  a  scalar  autocorrelation 
forms  a  min imura- delay  wavelet. 

(Kunetz  (1964)  has  proven  that  a  synthetic 
seismogram  which  includes  all  multiple  reflections  forms 
one  side  of  an  autocorrelation  function.  In  view  of  the 
results  obtained  above,  we  can  sharpen  his  result  to  say 
that  a  synthetic  seismogram  which  includes  all  multiples 
and  the  initiating  pulse  is  minimum  delay,  if,  and  only 
if,  the  initiating  pulse  is  minimum  delay  (Kunetz  took  this 
pulse  to  be  a  unit  spike,  which  is  certainly  minimum 
delay ) . ) 
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4.11  Factorization  Theorem 

An  autocorrelation  function  may  always  be  factored 
to  give  a  wavelet  a(z)  such  that  its  autocorrelation, 
a(z)  a#(z),  equals  rbe  original  autocorrelation.  In 
general  the  factorization  is  not  unique  but  it  may  be  made 
unique  by  requiring  that  a(z)  be  a  one-sided  invertible 
wavelet,  i.e.  minimum  delay.  Then  this  wavelet  is  the 
Robinson  canonical  minimum-delay  form  of  all  other  factori¬ 
zations.  These  properties  are  stated  more  rigorously  in 
the  factorization  theorem. 

Theorem  4.1-1  (Autocorrelation  Factorization) .  Let  r(z) 
be  a  real  scalar  autocorrelation  of  degree  +  m  .  Then 
there  exists  a  real  polynomial  (wavelet)  a(z)  of  degree 
m  such  that 

a)  r(z)  *  a(z)  a#(z)  . 

b)  a(z)  and  a-1(z)  an  both  analytic  inside 
the  unit  circle  Jz  |  <1,  i.e.  a(z)  is 
one-sided  invertible,  minimum  delay,  or  mini¬ 
mum  phase. 

c)  a(z)  is  unique  up  to  within  a  trivial  all¬ 
pass  system  multiplier,  i.e.,  if  b(z)  also 
satisfies  a)  and  b),  then  b(z)  =  p  a(z) 
where  p  is  a  constant  such  that  p  p«  1  , 
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d)  Any  factorization  of  the  form  r(z)  *  c(z)c#(z) 
in  which  c(z)  is  not  invertible  is  given  by 

c(z>  -  a(z)  p(z) 

p(z)  being  an  arbitrary  regular  all -pass 
system. 

Since  the  proofs  of  parts  c)  and  d)  are  very  similar  to 
that  for  matrix-valued  autocorrelations  we  will  defer  the 
proof  of  those  parts  until  the  next  section  (also  see 
Robinson,  1963.  p.  179).  The  proofs  of  j  'rts  a)  and  b) 
consist  of  showing  that  a  factorization  with  the  needed 
properties  exist.  We  will  state  three  factorizations  here 
but  will  defer  again  until  the  next  chapter  for  the  dis¬ 
cussion  of  approxiamte  factorizations  since  the  scalar 
methods  are  Just  special  cases  of  the  matrix- valued  tech¬ 
niques. 

4.12  Methods  of  Factorization 

4.121  Woldian  or  spectral  analysis 

As  pointed  out  at  the  beginning  of  section 
every  zero,  a^;  of  the  polynomial  r(z)  is  associated 
with  a  zero  l/ii  .  Thus  if  we  choose  the  m  zeros 
a ^  i  ■  1.  ...,  m  which  fall  outside  the  unit  circle  to 
form  the  polynomial,  a(z),  then  this  polynomial  will 
certainly  concur  with  parts  a)  and  b)  of  the  factori¬ 
zation  theorem. 


75 


4.122  ,  4.123 


4.122  Kolmogorov 

If  we  liave  the  square-gain  (that  is,  the  cosine 
transform),  r(e"lu!).  of  the  system,  then  the  wavelet, 
a(z),  is  given  by 

a(z)  *  £  a.z1 

1=0  1 

-  expLi-  f  “TS*5-  loe  '“(e1'")  <*»| ,  I  z  |  <  1 

L  -tt  e  -z  J 

(Robinson,  1963b  or  Karhunen,  1949) . 

4.123  Zero-phase 

The  zero-phase  factorization  is  also  based  upon 
the  cosine  spectrum,  however,  it  does  not  produce  a  wave¬ 
let  that  satisfies  part  b)  of  Theorem  4.1-1.  If  we 
desire  the  wavelet  to  be  two-sided  and  symmetrical  then 
we  need  only  take  the  square  root  of  the  spectrum 
jt(e"iu))  =  v/r(e"iu))  .  This  wavelet  has  zero  phase. 

The  spectral  and  Kolraorgorov  factorizations  are 
equivalent  (Robinson,  1954).  The  spectral  technique  is  not 
a  good  computational  method  because  of  the  well  known 
difficulties  in  finding  the  zeros  of  a  polynomial.  The 
Kolmorgorov  technique  becomes  approximate  in  computer 
applications  since  we  must  compute  some  continuous  func¬ 
tions  digitally.  It  has,  however,  been  successfully 
applied  to  factorization  problems  (Galbraith,  1963). 
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4.2  One-Dimensional  Matrix  Au t ocorr ela t ions 

The  matrix-valued  autocorrelation  function  is 
very  similar  to  the  scalar  function. 

Let  R(z)  be  an  n  x  n  quasipolynomial  auto¬ 
correlation  matrix  of  rank  r,  then 

a)  R(z)  is  reverse-synnetric,  i.e.  R(z)  »  R*(zj. 

b)  Rfe1^)  is  non-negative  definite,  i.e. 

#  . 

b  R(z)  b  >  0  for  every  n  vector  b  and 
every  value  of  z  on  the  unit  circle. 

c)  The  determinant  of  R(z),  d(z)  -  Jet  R(z) 
is  reveree-syumetrical  d(z)  »  d#(z)  . 

d)  The  Smlth-McMlllan  canonical  form  satisfies 

D(z)  «  D#(z)  . 

e)  The  real  frequency  zeros,  i.e.  the  zeros  on 
the  unit  circle,  of  the  diagonal  elements 
of  D(z)  (and  of  d(z))  are  of  even  multi¬ 
plicity. 

Prouf.  Statement  a)  is  obvious.  If  we  let  X(z)  re¬ 
present  an  arbitrary  finite  process,  then 

R*(*)  -  (X(z)  X.(z))# 

-  X..(z)  X.(z) 

-  X(z)  X#(z) 

”  R(z) 
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Statement  b)  follows  If  we  note  that  on  the  unit  circle 
the  determinate  of  R(z)  is 


|x{eitu)  X.(eil")  |  -  |x(elro)  |  |x  (e-iu>)| 

-  |x(elu<)  |  |  X  (eiu)| 

>  0 

unless  R(z)  is  null. 

Statement  c)  follows  directly  from  a).  For  statements 
d)  and  e)  we  let  R(z)  ■  C(z)  D(z)  F(z)  be  the  Smith- 
McMillan  canonical  form  of  R(z)  .  Now,  since  R(z)  = 
R#(z),  C (z )  D(z)  F(z)  »  F#(z)  D#(z)  C#(z)  .  But  D(z) 

and  D#(z)  are  both  canonical  to  the  same  matrix  R(z) 
and  therefore  by  the  Smith-McMillan  Theorem  must  be  the 
same.  Thus  every  diagonal  element  of  D(z)  is  reverse- 
symmetric,  and  consequently  any  zero  z^  is  accompanied 
by  a  zero  1/z^  ,  However,  if  z^  is  a  zero  of  R(z) 
then  it  must  also  have  been  a  zero  of  X(z)  .  Since  X(z) 
Is  real,  it  has  a  real  canonic  form  D^(z)  and 

D(z)  -  D1(z)  (D2)#  (z) 


But  since  is  real,  every  root  z^  must  be  accompanied 

by  its  complex  conjugate  z^  .  Therefore  for  every  root 
|  z^|  =  1  on  the  unit  circle  we  must  have  four  roots  z^, 
l/z^,  zlf  and  1/  z^  .  But  z^  ■  1/  z^,  z^  «  1/z^  if 
|z^|  *  1  .  Thus  all  roots  on  the  unit  circle  must  occur 


in  pairs. 


Q.E.D. 
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4.21  Factorization  Theorem 

A  matrix-valued  autocorrelation  may  always  be 
factored  into  the  product  of  a  wavelet  with  its  reverse- 
transpose.  This  factorization  is  made  unique  if  we  re¬ 
quire  that  the  wavelet  be  one-sided  and  invertible,  i.e. 
minimum  delay.  This  review  is  stated  more  concisely  in 
the  Matrix-valued  Factorization  Theorem. 

Theorem  4.2-1  (Matrix- Valued  Autocorrelation  Factorization) . 
Let  R(z)  be  a  real  n  x  n  quasipolynomial  autocorrela¬ 
tion  matrix  of  rank  r  .  Then  there  exists  a  real  n  x  r 
polynomial  matrix  A(z)  such  that 

a)  R(z)  ■  A(z)  A#(z) 

b)  A(z)  and  A  *(z),  its  left  inverse,  are 
both  analytic  inside  the  unit  circle.  If 
R(z)  is  full  rank  and  non-degenerate,  A(z) 
is  minimum  delay. 

c)  A(z)  is  unique  up  to  within  a  real-orthogonal 
matrix  multiplier  on  the  right  (a  trivial  all¬ 
pass  system),  i.e.,  if  Ax(z)  also  satisfies 
a)  and  b),  then  A1(z)  -  A(z)  T  where  T 
is  r  x  r,  constant  and  unitary,  T  T#  -  lr  . 

d)  Any  non-minimum  delay  factorization  of  the 
form  R(z)  .  C(z)  C#(z)  in  which  C<z)  is 
n  x  m,  m  >  r,  and  polynomial,  is  given  by 
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the  Robinson  canonical  form 

C(z)  -  A(z)  [  lr  |  Cr>m_r]  P(z)  . 

P(z)  being  an  arbitrary  rational  regular 
m  x  m  reverse-unitary  matrix  (that  is,  P(z) 
is  an  m  x  m  all-pass  system) . 

The  proof  to  this  important  theorem  is  divided 
into  two  parts •  First  we  prove  parts  c)  and  d).  Then 
parts  a)  and  b)  are  proven  in  the  next  section  by 
demonstrating  factorization  algorithms  which  produce  wave¬ 
lets  having  the  given  properties.  Four  such  algorithms  are 
known.  Two  produce  A(z)  by  analytical  manipulations  and 
two  give  A~^(z)  by  approximate  techniques. 

Proof.  Consider  statement  d)  first.  Let  C(z)«A(z)Q(z) 
where  A(z)  satisfies  a)  and  b).  Then 

C (z)  C#(z)  =  A(z)  Q(z)  Q#(z)  A#(z) 

=  A(z)  A#(z) 

Q(z)  Q#(z)  =  lr 

where  Q(z)  «  A’^z)  C(z)  is  obviously  analytic  inside  the 
unit  circle,  i.e.  P(z)  is  an  arbitrary  mx  m  reverse- 
unitary  matrix  that  incorporates  Q(z)  in  its  first  r 
rows;  i.e.. 


Q(z) 


•>  m 

1„  (L.  „  „  P(z)  . 

r  r,m-r  ' 

. 
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Now  let  us  consider  statement  c).  Let  A (2) 
and  A^fz)  be  two  matrices  satisfying  a)  and  b),  and 
let  A(z)  -  A1(z)  Q(z)  .  Then 

A(z)  A#(z)  -  A1(z)  (Aj).  (z) 

-  *!<*)  «(z)  Q.(z)  (A1),  (z) 

Q(z)  Q.(z)  -  lp 

where  Q(z)  «  A11(z)  A(z)  Is  analytic  Inside  the  unit 
circle.  But  we  also  have  Q(z)  -  (Aj),  (z)  A^z)  and  It 
is  therefore  analytic  outside  the  unit  circle.  By  Theorem 
3,3-6#  Q(z)  is  a  constant  real  orthogonal  matrix, 

4.22  Analytic  Factorization  Methods 

Both  of  the  analytic  factorizations  depend  upon 
the  factorization  of  an  elementary  autocorrelation  matrix. 
We  will  discuss  this  technique  first.  The  algorlttmi  was 
first  presented  by  Oona  and  Yasuura  (1954,  pp.  125-177) 
and  later  expanded  up^n  by  Youla  (1961,  pp.  176-178)  for 
paraconjugata-hermitian  matrices.  The  following  statement 
has  been  altered  to  account  for  the  properties  of  reverse- 
symmetric  matrices. 

4.221  Elementary  autocorrelation  matrix 

Consider  an  r  x  r  positive  elementary  quasi- 
polynomial  reverse-symmetri  matrix,  i.e.#  an  elementary 
autocorrelation  matrix,  R(z)  .  Because  of  t.,  positive 
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nature  of  R{eilt),  all  its  diagonal  elements  are  reverse- 
synmetric  and  positive  on  the  unit  circle.  Let 

£  q2  — *  •  ••*  <  Qr  be  the  maximum  degrees  of  the  diag¬ 
onal  entries  arranged  in  non-decreasing  order.  Since 
R(z)  is  reverse-synmetric,  the  q's  are  non-negative 
integers.  Again  invoking  the  positive  character  of 
R(eiU)),  it  follows  that  no  element  in  R(z)  has  degree 
exceeding  qp  .  Thus  qp  =  0  if  and  only  if  R(z)  is  a 
constant  symmetric  positive-definite  r  x  r  matrix,  in 
which  case  it  can  be  written  as  AA  by  a  number  of  stand¬ 
ard  techniques.  Excluding  this  relatively  trivial  situa¬ 
tion,  we  will  assume  qp  >  0  . 

We  begin  by  interchanging  the  rowB  and  columns 
of  R(z )  so  as  to  make  its  diagonal  elements 
(R)J1,  (R)22*  •••*  (R)rr  Possess  the  degrees  qi,  q2,  ..., 
qp,  respectively.  Call  the  rearranged  matrix  R^fz)  . 

Then  there  exists  a  permutation  matrix  K  such  that 

R2  (z)  =i  K  R(z)  K’  (4 .2-1) 

R^  is  also  elementary,  reverse-symmetric  and  positive. 

Next  we  force  each  diagonal  term  to  have  degree 
qr  .  Let  us  begin  by  defining  a  non-increasing  sequence 
of  non-negative  integers  <r^,  <7^,  ...,  <rv  by 

<r±  =  Qr  “  q^  i  =  1,  2,  ...,  r 

(4.2-2) 
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and  the  r  x  r  diagonal  matrix  H(z)  by 

r  k.  k  cr‘ 

H(z)  =  dlag  1^(1  -  z  *)  1,  (i  -  z  *)  ,  ...,  (1  -  z  r) 

(4.2-3) 

where  =  +  1  chosen  so  that  the  degree  of  the  non- 
diagonal  ternu*  do  not  exceed  qr  .  Note  that  0~T  *  0  . 
The  r  x  r  matrix 

Hg (2 )  -  H(z)  n1(z)  H,<«)  (4.2-4) 

is  quasipolynomial,  reverse-symmetric  and  positive. 
Moreover  all  of  its  diagonal  elements  have  the  same  degree 
qp  .  Since  R^  is  elementary,  it  is  clear  that 

Det  Rg (z)  *  O(z^)  (4.2-5) 

where 


^  -  ^1  +  ®2  +  •••  +  rr  (4.2-6) 

But  from  equation  4.2-2 

<r  <  (r  -  1)  qp  .  (4.2-7) 


Rg (z)  may  be  written  in  expanded  form  as 

/  %  “Qp 

Rg(z)  -  z  +  ...  +  I”  z  +  TQ  +  Tx  z  + 


...  +  T  z 
qr 

(4 .2-8) 


where  the  T’s  are  constant  r  x  r  matrices.  The  impor¬ 
tant  observation  is  that  T  is  singular,  i.e.  Det  T  ■ 

Qr  Qr 

0  for  otherwise  equation  4 ,2-8  would  yield 
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Det  R2(z)  =  G(z  r) 


(4 .2-9) 


which  contradicts  equations  4.2-5  and  4.2-7.  This  deduc¬ 
tion  implies  that  T  contains  a  principal  minor  Q  of 

4r 

order  s  x  s  which  is  non-singular  and  such  that  the  minor 

“3  created  by  adding  the  (s  +  l)th  row  and  column  to  G 

is  singular.  Thus  we  may  add  a  linear  combination  of  the 

first  s  rows  of  T  to  the  (s  +  l)th  row  and  the 

qr 

same  linear  combination  of  the  first  s  columns  to  the 

(8  +  l)th  column  such  that  (T  )  ^  is  reduced  to 

r 

zero  and  no  other  diagonal  term  is  affected.  Hence  for 
the  correct  choice  of  a  constant  r  x  r  non-singular 
matrix  K^, 


Cr  '  Kl\Kl 


(4.2-10) 


has  a  zero  element  in  the  (s+1,  s+1)  place.  From  4.2-8 


•%_  i 


R3(z)  =  Kx  R2(z)  K1=^  (Kx  T1  Kx)z 


(4.2-11) 


has  a  diagonal  element  in  the  (s+1,  s+1)  position  of 


degree  <  qp  . 


The  matrix 


R„(z)  =  H-1  (z)  R3(z)  H^z) 


.2-12) 


is  reverse- symmetric,  positive  and  elementary.  According 
to  the  definition  of  R0  (see  equation  4.2-4)  (R0).  _  is 
divisible  by  (1  -  z  ^)  (1  -  z  m)  m,  and  according 
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to  the  definition  of  R^(z)  (see  equation  4.2-11)  and  the 
definition  of  K^,  R^(z)  differs  from  Rg(z)  only  in  its 
(s  +  l)th  row  and  column.  More  specifically, 

^Vt.s+l  "  '^t.s+l  +  ^2^4,1  2’  r)> 

(4.2-13) 

the  c’s  being  scalars.  By  construction  6r^  >  <r^  >  > 

<7*1,  thus  every  term  on  the  right-hand  side  of  equation 

k_.  T  o-l.  !  -k  <7“ 

4.2-13  is  divisable  by  (1  -  z  8+1 )  8+1  (1  -  z  ra)  m, 

(k  =  1,  2,  ...,  r)  .  The  same  considerations  apply  to  the 

(s  +  l)th  row,  whence,  for  all  l  and  m,  (R~),  is 

k  o~  -k  <r~  3  m 

divisible  by  (1  -  z  *)  1  (1  -  z  ^*)  m,  and  R4(z)  Is 

a  quasipolynomial  matrix.  Since 

Det  Rjj(z)  *  DetfKj2  K2)  Det(R(z))  •  constant, 
R^(z)  is  elementary. 

But  R^(z)  is  simpler  than  R^z)  because  the 
degree  of  its  (o+l,  s+1)  entry  is  at  least  1  less  than 
the  same  entry  in  the  latter  matrix,  while  all  other  corres¬ 
ponding  diagonal  elements  have  the  same  degree  as  before. 
Consequently,  after  one  cycle  of  the  algorithm, 

R(z)  =  r,(z)  Rji  (z )  (C.)#(z) 

(4.2-14) 

where 

G1(z)  -  K"1  H"1^)  K^1  H(z) 
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1b  an  elementary  polynomial  matrix  and  R^(z)  1b  at 
least  1  deg:  less  than  R(z)  . 

We  now  replace  R(z)  by  R^(z)  and  repeat  the 
algorithm.  After  a  maximum  of  rqr  cycles  R(z)  is  re¬ 
duced  to  a  constant  symmetric  positive-definite  matrix 
Rij  *»  CC*  t  so  that  finally 

R(z)  -  A (z)  A,(z) 

where 

A (s)  =  G1(z)  G2(z)  G  (z)  C  . 

This  factorization  does  not  guarantee  that  A(z) 
is  one-sided.  This  is  because  of  the  ambiguity  in  the  defi 
nition  of  H(z)  .  Tc  the  author's  knowledge  ao  one-sided 
factorizations  exist  for  cases  in  which  this  algorithm 
does  not  give  a  right-sided  factorization.  For  example, 
the  elementary  autocorrelation 
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Figure  4.2  -  Is  Elementary  autocorrelation  matrix 
factorization. 
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but  no  right-sided  form  has  been  found. 

The  important  steps  in  this  reduction  are  shown 
in  Figure  4.2-1. 

Example  4.2-1.  Let  us  consider  the  elementary  autocorrelar 
tion  matrix 

R(z)  *  Uz-1*  6-  2z  -4z~2+l4z~1-l4+  4z  1 

L  4z“1-l4+l4z-4z2  8z“2-32z“1+50-32z+8z2J 

V/e  will  follow  the  steps  of  this  factorization  in  detail. 
Recursion  1. 

Since  the  degrees  of  the  diagonal  terms  of  R(z) 
are  already  in  ascending  order  we  may  skip  the  first  step. 

Thus 

R-j^Cz)  =  R(z) 

Next  we  make  all  of  the  terms  have  the  same  degree  by  form¬ 
ing  the  product 

rJ1^(z)  =  R^1)  =  diag[l-z,  l] 

-  “ 

=  2z“2-10z"1+l6-10z  +2z2  -4z“2+18z"1-28+18z-4z2 

-4z2  +18z“1-28+18z"1-4z2  8z“2-32z_1+50-32z+8z2 

-  -J 

2 

The  matrix  for  the  z  terms  is 

T21)  *  ^  2  “4J 

| 
| 
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Thus  the  diagonal  term  (T2)2  2  can  be  re<*uced  to  zero 
by  the  product 


R^(z)  .  rJ1*  (K^1^ ) '  where 


[;  3 


=  -2Z*1  -  lOz”1  +  16  -  lOz  +  2z2  -2z-1  +  4  -  2z 


-  2z“1  +  4  -  2z 


Finally  we  remove  the  H  multipliers  to  obtain: 


R^(z)  -  1  Ri1) 


-2z-1  +  2 


•£z  ^  +  6  -  2z 
2  -  2z 


These  steps  can  now  be  combined  so  that 


R(z)  «  G1  R^1'  (Gx) 


where 


C,  -  1  lei1*  1 


-2+2 z  1 


Recursion  2, 


Our  beginning  point  Is  the  matrix  RjJ1'  from  the 
last  recursion.  This  time  we  oust  exchange  the- positions 
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of  the  diagonal  terms  so  that  they  will  have  ascending 
degrees; 


r(2)  a  K(l)  R^l)  K(l) 1 

where  K  =  ^0  1^ 

-  f  2  2  -  2z 

L-2Z"1  +  2  ~2z_1  +  6  -  2z 

Here  we  must  now  multiply  by  H^2^  =  diag  [(l-2_1),l]  in 
order  for  the  off  diagonal  terras  to  be  the  same  degree  as 
the  diagonal  terras. 

rJ2)(z)  *  R^2^  Hp) 


-2z'2  +  4  -  2z 

-2z‘2  +  4  -  2z 

-2z_1  +  4  -  2z 

-2z“*  +  6  -  2z 

Prom  the  coefficients  of  z  we  may  select  k|2^  such  that 
the  degree  of  (r|2^)2  2  is  reduced, 

r|2'  =  k|2>  R^2>  where  K^2*  =  ^  1 

-2Z'1  +  4  -  2z~*  0  1 

0  2  J 

and,  proceeding  as  before,  we  have 

Rjj2^  «  (H^)”1  R^  (h£2^  )_1 

2  Ol 

0  2 
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This  Bet  of  operations  may  also  be  grouped 


C2  »  K*2)  1  1  JC-j2)  1  H^2^  . 

Now  we  see  that  the  factorization  of  R(z)  ia  given  by 
R(z)  -  A(z)  A#(z) 

where  A(z)  -  r2 


ST 


■  z"1  +  1 
2Z*1  -  3  +  2z 


1 

-2  +  2z 


4.222  Spectral  analysis 

We  will  begin  by  illustrating  the  decomposition 
for  an  n  x  n  full  rank  (r  -  n),  non-degenerate  (p-2mn 
zeros  In  the  determinant,  where  m  is  the  greatest  number 
of  zeros  in  any  of  the  quasipolynomial  elements)  auto¬ 
correlation. 

Let  us  assume  that  statement  a)  is  true.  We 
begin  then  by  examining  the  latent  zeros  and  vectors  of 
R(z)  in  terms  of  those  of  A(z)  .  The  latent  roots  as 
specified  by  the  determinat 

|R(z)|  .  |  A(z)  A.(z)  | 

-  |A(z)  |  |a(z_1) |  .  (4.2-15) 

are  zy  1  ■  1,  2,  ....  p  where  the  z,  are  roots  of 
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A(z)  .  The  latent  vectors  are  either 

Adj  R(z1)  =  Adj  (A (Zj^)  A*(zi) ) 

-  Adj  A#(zi)  Adj  A(zi) 

*  Adj  A*^)  (4.2-16) 

=  2Li  v[ 
or 

Adj  Rfz^1)  *  v1  ui  Adj  A(zi*) 


If  we  choose  the  p  zeros  outside  the  unit  circle  we  will 
satisfy  condition  b).  These  zeros  and  their  associated 
vectors  v,  may  be  used  to  construct  7T (z )  according  to 
either  of  the  two  methods  illustrated  in  the  proofs  to 
Spectral  Theorem  or  the  Spectral  Corollary  (Section  3.32). 
We  must  now  determine  the  constant  multiplier  AQ  from 
the  autocorrelation 


R(z)  «  E(z)  AQ(A0)#  £#(z) 
AqAq  *  ®”1(z)  R(z)  (z)  . 


(4.2-17) 


Thus  we  can  only  determine  AQ  to  within  a  real  orthogonal 
multiplier. 


The  factorization  used  above  is  similar  in  intent 
to  the  Woldian  factorization  for  scalar  autocorrelations. 
The  °quence  of  operations  Is  illustrated  In  Figures  4.2-2 
or  4.2-3. 
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t  i 


M  «  -  I 


r»gnd«d  down 


tt  m«  Figure  4#2"1 


Figure  4.2  -  3*  Spectral  factorization  of  a  matrix  auto* 
correlation  according  to  Corallary  3.3-1. 
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If  R(z)  is  full-rank  but  degenerate,  p  <  nm, 
then  the  factorization  is  not  complete  when  we  reach 
equation  4.2-17.  In  this  case  we  will  have 

S(^)  *  Tr1^)  r(z)  ^(z) 

where  I? ( z)  is  an  elementary  autocorrelation  matrix  which 
must  be  factored  according  to  the  method  of  Section  4.221 
to  give 

K(z)  -  B0(z)(B0).(z)  . 

The  complete  factorization  Is 

A(z)  «  B(z)  <?0(z)  . 


If  R(z)  is  not  full  rank,  then  the  factorization  must  be 
done  in  terms  of  full  rank  submatrices  of  R(z)  .  Thus, 
we  partition  R(z)  symmetrically  about  the  main  diagonal 
such  that  each  r1  x  submatrix  R^Cz)  is  full  rank. 
For  example 


•5 


Each  of  the  diagonal  subma trice  (R)^,  is  then  factored 
according  to  the  spectral  theorem  technique  (given  earlier 

A 

in  this  section)  to  obtain  the  r^  x  r^  matrices  G^fz)  • 
Now  form  the  matrices 


G(z) 


A 

1 

coo 

A 

Gn 

0-  „  v. 

ri’  r  ^1 

c22 

°r2,  r-r2 

* 

(4.2-19) 


and  the  left  inverse 


H(z) 


C 


- 

*-l 

C11 

C22 

“-I 
c  A 

UAA 

uu 

l 

0_  _ 

p_rl'  rl 

o 

i 

-i 

ro 

ro 

1 _ 

(4.2-20) 


where  C  is  a  constant  diagonal  r  x  r  matrix. 
Now,  the  matrix 


H(z)  *  g“*(z)  r(z)  n*1( z) 


is  an  r  x  r  elementary  quasipolynomial  matrix  (this  will 
not  be  proven  here)  which  may  be  factored  according  to 
Section  4.221. 

Thus 


H(z,  =  CQ(z)  (Gq)#(z) 
A (z )  «  G  (z)  r0(z) 


i 


| 

I 


and 
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Example  4.2-2.  (Pull  rar'',  non-degenerate  case) 


Consider  the  autocorrelation  matrix 

R(z)  =  -2s”1  +  6  -  2z  -z_1  +  1 

1  -  z  -z"1  +  2  - 

_ 

We  begin  by  finding  the  latent  roots  z^  and  latent 
vectors  v^  and  w^  by  the  technique  outlined  in  the 
section  3*32. 

R(z)  -  2z“2  -  9Z*1  +  14  -  9z  +  2z2 

-  (1  -  2z) (1  -  z)(l  -  z“1)(l  -  2Z-1) 

Substituting  these  roots  into  AdJ  R(z)  we  find 


zero 

vector  w 

vector 

1/2 

(i.  i) 

(-1.  2) 

1 

(o,  1) 

(0,  1) 

1 

{0,  1) 

(0,  1) 

2 

(-1.  2) 

(1.  1) 

We  choose  one  of  the  roots  on  the  unit  circle  and  the 
root  outside  the  unit  circle  to  find  A  •  Using  the  no¬ 
tation  of  section  3 #32,  we  have 
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Thus  ?f(z)  =  I  -  V  z 

=  fl  -  1/2  z  1/2  z' 
0  1  "  z. 


and 


AqA£  =  ^_1(z)  R  (z )  ^(z) 


1 

1-z 

l/2z 

-2z”^+6-2z  -z-1+1 

(l-3/2z+l/2z2)2 

0 

l-l/2z. 

1-  z  -z_1+2-z_ 

5  1 
1  1 


-z^+l  0 

l/Sz"1  -l/2z_1+l 


o' 

1 


Consequent  ,y  we  find  the  desired  minimum  delay  wavelet 


A (z)  ^  tf(z)  Aq 


2-z  1 

0  1-z 


9B 


1 
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4,223  Smith-McMillan 

This  factorization  is  based  upon  the  Smith- 
McMillan  canonical  form.  The  algorithm  is  very  similar  to 
a  factorization  technique  for  paraconjugate-hermitian  ma¬ 
trices  given  by  You la  (1961).  The  system  is  quite  elegant 
in  its  conception  since  the  algorithm  is  independent  of 
the  rank  and  degeneracy  of  the  autocorrelation;  however,  a 
flaw  remai  so  that  for  some  cases  a  one-sided  factoriza¬ 
tion  cannot  be  guaranteed  even  for  the  full-rank  non¬ 
degenerate  matrix  autocorrelation. 

Before  doing  the  factorization  we  must  investi¬ 
gate  several  mor"  properties  or  quasipolynomial  matrices. 

Definition  4.2-11.  Let  Q(z)  be  an  n  x  m  rational  matrix 
of  normal  rank  r  ,  A  decomposition  of  the  form 

G(z)  =  P(z)  £(z)  Q(z) 

is  said  to  be  an  inner- standard  factorization  if 

a)  4(z)  is  r  x  r,  canonic  and  analytic  to¬ 
gether  with  its  inverse  in  the  entire  z 
plane  with  the  possible  exception  of  a  finite 
number  of  points  on  the  unit  circle. 

b)  P(z)  is  n  x  r  and  analytic  together  with 
its  left  inverse  Inside  and  on  the  unit  circle. 

c)  Q(z)  Is  rxm  and  analytic  together  with  its 
right  inverse  outside  and  on  the  unit  circle. 
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Interchanging  the  roles  of  P  and  Q  gives  rise  to  an 
outer- s tandard  factorization.  Obviously  any  inner-standard 


factorization  of  G(z)  generates  an  outer-standard  factori¬ 
zation  of  G’(z),  G~^(z)  and  G(l/z)  ,  For  example  G' (z)  = 

Q*  (z)  A(z)  P' (z),  etc. 

It  follows  from  the  Smith- Me Mil lan  Theorem  that 
any  rational  matrix  G(z)  possesses  an  inner-  and  outer- 
standard  factorization.  For,  let  G(z)  =  C(z)  D(z)  F(z) 
where  C  and  F  are  elementary  and  D  is  canonic.  By 
factoring  the  cp's  and  4’s  (see  Smith-McMillan  Theorem, 
section  3.3*0  appearing  in  the  diagonal  elements  of  D(z) 
into  the  product  of  three  quasipolynomials,  the  first 
without  zeros  j  z  |  <  1  ,  the  second  without  zeros 
|z|  /  1,  and  the  third  without  zeros  in  |z|  >  1,  it  is 
possible  then  to  write  D(z)  »  D+(z)  A(z)  D“(z):  D+(z) 
and  Its  inverse  are  analytic  in  |z|  <  1,  A(z)  and  A^(z) 
In  |  z  |  /  1,  and  D~(z)  and  its  inverse  in  |z|  >  1  . 

Now,  choosing  P(z)  =  C (z )  D+(z)  and  Q(z)  =  D~(z)  F(z) 
we  have  the  desired  breakdown. 


Lemma  4.2-1. 

Let  G(z)  possess  two 

right-stanaard 

factori- 

zations  G  * 

P1  ^1  ^1  *  Then* 

a) 

A(z)  =  Aj(z) 

b) 

P1(z)  =  P(z)  M*1(z) 

and  Q^(z)  = 

N(z)Q(z), 

where  M(z)  and  N~^(z)  are  any  two  r  x  r 


elementary  qua si polynomial  matrices  which 


transform  A(z)  into  itself ,  vis, 

M(z)  A(z)  N_1(z)  -  A(z)  . 

Proofs  We  have 

0  -  PAq  -  PjA^  (4,2-21) 

Then 

A"1  P"1  PA  -  Qx  Q"1  (4.2-22) 

By  definition  the  right  hand  side  of  equation  4.2-22  is 
analytic  in  |  z  |  >  1  .  Thus  Q^"1  is  analytic  in  the 
entire  z  plane.  According  to  equation  4.2-21  the  inverse 
of  QjQ”1  is  A 1  P  1  Pj  A^  *  QQ^1  and  is  therefore 
also  analytic  in  the  entire  z  plane.  By  Theorem  3.3-4 
Q1Q~1  is  therefore  an  elementary  r  x  r  quasipolynomial 
matrix  N(z)  .  Similarly  p"1!*  is  an  r  x  r  elementary 
quasipolynomial  matrix  M(z)  .  Prom  equation  4.2-21 

M(z)A(z)  N-1(z)  *  Ax(z)  . 

Since  A(z)  and  A^(z)  are  both  canonic,  A( *>)  ■  A^(z) 
by  the  Smith-McMillan  Theorem,  Thus 

M(z)  -  ^(z)  N(z)  SX{z) 

Qx(z)  *  N(z)  Q(z) 

Px(z)  =  P(z)  A{z)  lT^z)  A1(p) 

=  P(z)  M"1(z)  Q.E.D. 
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Corollary  4 .2-1.  The  canonic  matrix  A(z)  appearing  in 
either  an  inner-standard  or  outer-standard  factorization  of 
an  n  x  m  matrix  G(z)  of  rank  r(G)  is  equal  to  the 
r  x  r  identity  matrix  lr  if  and  only  if  G ( z )  is 
analytic  and  r{G)  is  constant  on  the  unit  circle.  In 
this  case,  if  P  Q  and  Q1  are  any  two  standard 
factorizations  of  G,  P1(a)  =  P(z)  N_1(z)  and  Q1(z)  = 

N(z)  Q(z),  N(z)  being  an  arbitrary  r  x  r  elementary 
quasi polynomial  matrix. 

Proof.  The  if  part  is  immediate.  The  analyticity  of  G(z) 
on  the  unit  circle  implies  that  all  of  the  denominator  terms 
of  A( z)  are  unity.  This  in  turn  leads  to  the  conclusion 
that  r(G)  is  constant  on  the  unit  circle  only  if  the 
numerator  quasipolynomials  in  A(z)  are  unity.  Thus 
A(z)  ■  lr  •  The  remaining  statements  are  consequences  of 
Lemma  4.2-1.  Q.E.D. 

Corollary  4.2-la .  If  G(z)  is  reverse-symmetric  then 

N(z)  =  M#(z) 

where  M(z)  is  any  r  x  r  elementary  quasipolynomial 
matrix  satisfying  A(z)  M#(z)  =  M(z)  A#(z)  . 

Proof .  Since  G(z)  =  G#(z),  Q*(z)A(z)  P#(z)  is  also  a 
right  standard  factorization  of  G(z)  by  arguments  similar 
to  those  used  for  theorem  3*3“4.  Thus,  according  to 
Lemma  4.2-1 


±02 


P*(z)  =  N(z)  Q{z) 

Q„(z)  -  P(z)  M*1(z) 

P.(z}  -  N(z)  M;a(z)  P,(z) 

Since  P#(z)  has  a  right  inverse, 

N(z)  -  M#(z) 

and  according  to  Lemma  4.2-1 

Zi(z)  M#(z)  *  M(z)  4#(z)  .  Q.E.D. 

The  factorization  algorithm  discussed  here  is 
based  upon  the  Smith- Me Mil lan  canonical  form  for  the  auto¬ 
correlation  matrix.  Unfortunately,  because  of  the  arbi¬ 
trariness  of  the  sequence  of  steps  involved  in  finding  a 
particular  realization  of  the  Smlth-McMlllan  canonical 
form,  the  solution  Is  not  unique.  The  solution  matrix 
X(z)  Is  not  one-sided  (and  therefore  not  analytic  Inside 
the  unit  circle).  This  matrix  X(z)  will  differ  from  the 
proper  answer  by  a  unitary  matrix. 

Step  1.  Reduce  the  matrix  R(z)  to  its  Smith- Me  Mil lan 
canonic  form.  Since  R(z)  Is  a  quasipolynomial  matrix, 
this  procedure  Is  a  standard  but  arbitrary  one  as  Illus¬ 
trated  in  section  3*34.  Thus  we  will  have 

R(z)  -  C(z)  D(z)  F(z)  . 

Step  2.  According  to  Theorem  4.2-1,  D(z)  is  of  the  form 
that  it  may  be  factored  as 
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D(z)  =.  D+(z)  A(z)  D- (z) 

-  D+(z)  A(z)  Dt(z) 

where 

(1)  D+(z)  is  r  x  r,  diagonal  and  analytic, 
together  with  its  inverse  (D+)-1(z)  for 

z|  <  1  . 

(2)  A»(z)  =  A(z)  =  e(z)  Q#(z)  in  which  all 
diagonal  elements  of  0(z)  are  reverse- 
symmetric.  Furthermore,  A(z)  is  canonic 
and  non-zero  for  jzj  /  1  . 

Let 

P(z)  =  C(z)  D+(z) 

Q(z)  =  Dt(z)  F(z) 

Then  we  have  an  inner-standard  factorization 
R(z)  =  P(z)  A( z }  Q(z) 

Step  3*  Now  we  wish  to  factor  A(z)  •  Since  R(z)  is 
reverse-symmetric,  a  second  left  standard  factorization  is 

R(z)  =  Q#(z )  A*(z)  P#(z) 
and  according  to  Lemma  4,2-1  and  its  Corollaries 

Q*(z)  =  P(z)  M”1 (z)  (4.2-23) 

where  M”^(z)  is  an  r  x  r  elementary  quasipolynomial 
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matrix  such  that 

^(z)'1  M(z)  A(z)  »  N(z)  (4.2-24) 

Is  also  qua si polynomial , 


Thus  we  may  write 

R(z)  -  P(z)  M‘1(z)A(z)  P.(z) 

-  pee'1  m'1  e  e#  p# 

or 

9 _1  p"1  R  p;1  9;1  .  9"1  m”1  9 

Hence 


(4.2-26) 


(4.2-26) 


f5(z)  a  9"1(z)  M-1(z)  9#(z) 


is  r  x  r,  reverse-symmetric  and  non-negative  on  the  unit 
circle  (by  the  properties  of  equation  4.2-26).  Actually 
we  can  say  a  good  deal  more.  Let  us  write  equation  4.2-24 
in  terms  of  its  elements: 


(M)rk(z) 


CgjfcU) 

(A)rr(z) 


(H)rlc(z) 


(®)pr(z) 


Since  each  element  must  be  quasipolynomial 


(M)rk(z) 


(9WZ> 

(9)rr(z) 


must  also  be  quasipolynomial.  Thus  7f(z)  is  a  quasi- 
polynomial  matrix.  But  $(z)  »  M-1(z)  *  constant,  i.e,, 


7 


105 


Wil^Ui]llliilKI^W1NllMWllWHilllWWH^lll|im,,www,:MtHU,>,>>i:>:‘*H>>K,H<'',',,',,>H .  . . . . . — . . . . * 


$(z)  is  a  positive  reverse-symmetric  r  x  r  elementary 
quasipolynomial  matrix.  In  Section  4.221  we  demonstrated 
that  such  a  matrix  is  factorable  as 

$(z)  =  S (z )  S#(z)  , 

S(z)  being  an  r  x  r  elementary  quasipolynomial  matrix. 
After  this  is  achieved,  a  factorization  for  R(z)  is 
obtained  as  R(z)  =  7C(z)  A#(z)  with 

X(z)  =  P(z)  P(z)  S(z) 

-  C(z)  D+(z)  f (z)  S(z)  . 

where  A(z)  differs  from  the  desired  factorization  A(z) 
by  a  unitary  matrix.  By  straight  forward  algebra 

X(z)  £#(z)  =  C  D+  0  S  S#  B*  Dl  C* 

=  C  D+  1  p*  D*  C# 

=  P  M*1  A  P# 

=  Q.AP, 

=  R 

The  pertinent  computational  steps  involved  in 
this  algorithm  are  illustrated  in  Figure  4.2-4. 

The  advantage  of  this  factorization  is  that  the 
degenerate  and  singular  autocorrelation  matrices  need  not 
be  treated  as  special  cases  (as  contrasted  to  the  spectral 
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L  R(z>  •  R^ i ) 


RoMkco  to  Smith  - 
McMillan  canonicol 

f  Of  1*1 


Cfz)DIZ)r(z)  >R(zi 


OMam  xwar-otawMarM 

factorization 


Rtf)fti)Vi>0(i)*R(i) 


CoRtRwM  rtta  olxmootory 
matrix 

M'1  •  R*'0, 


.(If:) 


Comfuto  Mia  oiomaitiary 
ow  toe  or  rotation  matrix 

M(i)«r'M*'i 


pi)l(l)Mll)^(i)||t)iR(|) 


All )  Aa(z ) 


Factor  Mm  otomoatary 
aw  toco?  ratal  km  matrix  M 

ss.  •  5 


AUTOCORRELATION 
FACTORIZATION 
A ( z)  •  P(z)flz)Slz) 


MSS,«,PW  *R(z) 


**  See  Pig*  4,2-1, 


Plt^ure  4,2  -  4:  Smith- McMillan  factorization  of  a 
matrix  autocorrelation. 
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approach) .  If  an  algorithm  can  be  found  for  forming  the 
Smith-McMlllan  canonical  form  so  that  the  factorizations 
must  be  one-sided  then  this  formulation  would  become  more 
important  than  the  spectral  approach. 


Example  4,2-3.  As  in  Example  4.2-2  we  will  consider  the 
autocorrelation  matrix 


R(z)  = 


-2z  ^  +  6  -  2z  -z^+1 

1  -  z  -z'1  +  2  - 


£ 


Step  1  of  the  factorization  consists  of  reducing  R(z)  to 
the  Smith-McMlllan  canonical  form.  Since  this  process 
was  Illustrated  in  Example  3 • 3“3  we  will  merely  give  a 
particular  result: 


R(z) 

where  C(z) 


D(z) 


and  F(z) 


C(z)  D(z)  F(z) 

r  n 

1  o 

2  -  9/2z  +  7/2z2  -  z3  1_ 

diag  [l,  z"2  -  9/2z_1  +  T  -  9/2  z  +  z2] 

[-2z_1  +  6  -  2z  -z"1  h  1  “ 

4z  -  2z2  z 


Notice  that  the  factorization  has  been  made  in  such  a  way 
that  C(z)  is  one-sided. 


Step  2  consists  of  forming  the  left-standard  factorization 
from  the  Smith-McMlllan  canonical  form.  We  write 
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D(z)  =  dlag  [l,  (1  -  l/2z)(l  -  z)  2  (1  -  z"1)(l  -  1/pz"1)] 
=  D+  9  0#  D* 

where  D+  =  diag  [l,  (1  -  l/2z)J 

0  =  v'T  diag  [l,  (1  -  z)]  . 

Now,  the  left  standard  factorization  R  =*  P  A  Q  is  given 
by  setting 

P  »  C  D+ 

.[’  2  “  1 

[2  -  9/2z  +  7/2z  -  zJ  1  -  l/2zj  , 

Q  -  I>t  P 

[-2z_1  +  6  *  2z  -z"1  +  ll 

-2  +  5z  -  2z2  -1/2  +  zj 

and  A  »  0  0  # 

Step  3  involves  extracting  the  elementary  reverse-symmetric 
polynomial  'H(z)  from  the  left-standard  factorization. 

We  have 

R  -  P  A  Q 

«  P  M*1  A  P# 

where  m’1  -  P_1  Q#  (see  equations  4.2-23  and  4.2-25). 

Thus  we  compute  M’1 

,  f-2z-1+  6-  2z  -2z*2+  5z-1-  2 

*  a 

L  8z_1-  9fl4z-  9z2+8z3  4z-2-16z"1+25-16z+4z2. 

Now,  we  also  had 
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4.223 


r  =  pee"1  m"1  e  e# 


=  p  e  e#  p# 


where 

H  = 

e"1  m-1  e 

r  -1 

-2  „  -1 

-2z  •>  6-  2z 

-2z  +  7z  -  7+  2z 

„  2z_1-  7+  7z-2z2 

4z  2-l6z  1+25-16z+4z 

This  is  indeed  an  elementary  reverse-symmetric  polynomial 
matrix.  The  next  3tep  is  the  factorization  of  this  matrix 
into  the  form  $  =  S  S#  .  A  very  similar  elementary  auto¬ 
correlation  matrix  was  factored  in  the  Idxample  4.2-1.  The 
result  is 


S(z) 


i  1 

-2+2z  J  . 


The  factorization 

R(z)  = 
vjhere  A(z)  = 


of  R (z )  is  given  by 
A(z)  A#(z) 

P(z)  P(z)  S(z) 

1 

1/2  -  l/2z 


1 

l/2z  -  l/2z2 


This  is  not  a  one-sided  factorization  as  we  had 
obtained  in  the  spectral  decomposition  example;  however* 
this  solution  can  be  forced  to  be  one-sided  by  post- 
multiplying  it  be  the  proper*  unitary  matrix. 
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4.23  Approximate  Factorization  Methods 

The  factorization  methods  outlined  in  the  last 
two  sectioiBare  exact,  but  are  difficult  computationally. 

The  spectral  apprrach  suffers  from  the  well-known 
difficulties  of  determining  the  zeros  of  polynomials. 

The  Umith-McMillan  canonical  form  approach  is  complicated 
and  has  not  yet  been  refined  to  give  one-sided  factorizations. 

In  this  section  we  will  discuss  two  approximate 

schemes  for  determining  the  Taylor  expansion  of  the  inverse 
'x-l 

operator  a  (z)  from  the  autocorrelation  R(z).  Both  of 
these  techniques  depend  upon  the  fact  that 

7I'1(z)  R(z)  .  A0A0#E#(z) 

is  one-sided  (specifically,  right-sided).  Thus  if  we  have 
an  approximation 

i<1(z)  sa  rr1 

then  we  can  improve  the  approximation  be  examining  the 
non-zero  right  side  of 

4i(z)  R(z)  *  l(z)  . 

The  first  technique  is  a  recursive  method  that 
may  be  associated  with  least-squares.  It  was  advanced 
independently  by  Robinson  (1963&)  and  J,  P,  Burg  (personal 
communication)  based  oi  the  work  of  Levinson  (1947). 


Ill 
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The  second  technique  is  an  iterative  method  based 
on  the  vector  projections  of  linear  algebra.  It  was 
developed  by  Wiener  and  Masanl  (1957  and  19'  5  and  by 

Masanl  (i960).  Wunsch  (1965)  has  also  published  a 
heuristic  interpretation  of  the  projection  technique. 

Both  computational  schemes  have  been  programmed  and  tested 
for  computational  efficiency.  For  all  cases  tested,  the 
projection  technique  was  an  order  of  magnitude  slower  than 
the  least-squares  recursive  method. 


4.231  Least-squares 


The  approximate  least-squares  wavelet,  A„  , 

M 


of  degree  M 


Al\  =  Vm  +  *1,M  2  +  *•*  +  ^M, M  2 


M 


has  the  properties  that 


a)  Vm  =  *n  '  and 

b)  £ ,  M  =  0  is  1,  • .  • ,  M  where 

^jyj(Z  )  =  ^M(z)  r(z)  • 


M’ 


(4.2-27) 


If  we  write  out  the  equations  for  ^  i  =  1,  ..., 


M 


^M,M  R-r4fl  +  *•*  +  ^1,M  R0 

A  M,  M  R-Mf  2  +  * • *  +  A 1 , M  R1 


=  ’mo,m  ri 
“  "*0,M  R2 


(4.2-27a) 


R0  +  *  **  +  ^1,M  RM-1 


"  ”40,M  hm 


we  see  that  this  defines  a  set  of  nM  simultaneous 
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equations  which  can  always  be  solved  for  ^M(z)  •  We  will 
provide  a  recursive  technique  for  extending  the  length 
(degree)  of  ><M(z)  without  directly  resolving  the  set  of 
simultaneous  equations  given  above. 


For  the  recursion  we  sill  need  a  second  wavelet 

*  B0,Vl  +  *1,M  z  +  +  zM 


that  has  the  properties 


a>  a0,M  "  Xn  *  and 
b>  ?1.M  -  0  1  ’  -1 

,  ...»  -M  where 

*M(l/z)  R(z)  » 

*M<Z) 

(4.2-28) 

and  where 

^i,M  1  *  “1' 

-M  is  given  by 

*1,M  R-Mfl 

+  •••  +  *M,M  ro 

“  ~bo,m  r-m 

*1,M  R-Mf2 

• 

• 

+  **•  +  *M,M  R1 

• 

• 

"  “*o,m  r-hh 

•  • 

•  • 

(4.2-28a) 

• 

^1,M  R0 

• 

+  ...  +  M  RM-1 

•  • 

“  "5o,m  r-i 

In  the  spectral  factorization  of  the  autocorrelation 
(see  section  4.222)  we  obtained  the  minimum  delay  wavelet 
£(z)  by  choosing  all  of  the  roots  outside  the  unit  circle. 

We  could  also  have  formed  a  maximum  delay  wavelet  B(z) 
by  choosing  all  of  the  roots  inside  the  unit  circle  such 
that  R(z)  =«  B(z)  B0  Bq  B#(z)  .  Thus  B#(z)  is  minimum 
delay.  The  wavelets  ^M(z)  and  ^M(z)  are  the  least- 
squares  approximations  to  the  wavelets  ?T*(z)  and  S^(z). 


I 

-i 
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It  can  be  shown  (Robinson,  1963b)  that  4^(z)  and  i?M(z) 
are  also  minimum  delay. 


Notice  that  if  weight  and  add  ^(z)  and  £^(z)> 

we  find 

K  >»M(z)  +  z1*1  15  5„(  1/z)  R(z)  =  K  CM(z)  +  z^1  15  ?„(z)  . 

Thus,  if  we  choose  K  and  R  so  that 

a)  K  =  ln  ,  ^  ^0,M  =  0 


or 


define  K&^M  *  ^0,M  9 


+  l 


r  1 

'0,1 

-  0 


b)  K  =  ln  ,  K  C0^M  . 

define  K^M  =  C0,M  ' 

we  find  a  recurrence  relationship 


+  zltl  Ka,«  S«^z* 


Wz>  =  Vz>  +  zKfl  “b.M  V1^) 


These  polynomials  are  multivalued  counterparts 
of  the  polynomials  orthogonal  of  the  unit  circle  treated 
by  Geronlmus  (I960)  and  Szego  (1959) • 


Likewise  note  that 


Mfl^ 

zMfl  Ka,M 

-  Vz>  + 

z  CM^ 

There  are  two  other  relationships  that  are 
important  computationally.  First,  if  R(z)  is  symmetric 
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(see  Section  4.31)  then  ^M(z)  *  ^M(z)  ^or  obvious  reasons. 
Second,  for  all  cases 

CWfl,M  *  ^-M-1,M  * 


Proof.  (According  to  J.  P.  Burg,  personal  communication.) 
We  first  map  equations  4.2-27a  and  4.2-28a  into  matrix 
notation: 


[^M*  AVl,W  °] 

[°'  Bt\,n'  so,m]  bmh 

where 


[e0,M*  0>  °*  CHH,m] 

P-M-1,M*  °'  - -  °'  *0,M] 

(4.2-29) 


RMfl 


•  •  • 


RWfl 


The  solution  to  the  next  recursion  will  then  give 


.eo,w-i’  °*  •••’  °] 

*M,Mtl'  *o,w-i 

“wn  * 

0,  0,  *0>Mfl] 

(4.2-30) 

Equations  4.2-30  show  that  the  first  and  last  rows 
of  are 


C0,Mfl  [*0,MH'  *Mfl,Mfl] 

^0,Mtl  *  *  *  *  ^0,  Mfl]  • 


lib 


i 

I 


i 
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Since  R  is  symmetrical,  R  ^  is  also  symmetrical,  and 
therefore 

(e 

(note  that  ^  and  ^  must  be  non-singular  for 
all  M  since  R  is  non-singular.) 


-1 

0,Mfl 


Mfl 


V  _  r1 

,Mfl )  ~  C0, 


Mfl  *MH,Mfl 


(4.2-31) 


if  RMH 
^Mfl 


There  exists  an  n  x  n 
=  Q  ,  then  =  0  . 
.  However,  since 


matrix  Q  such  that 
If  RMfl  =  Q  * 


B  -  Tf  1 1  •  “ 

"  60,Mfl  O.Mfl 

(from  equation  4.2-31)  if  =  0,  then 

and  ^Mfl^z)  =  (z)  .  For  an  arbitrary  we  can 

write 


R 


Mfl 


-  <R«fl  -  «)  +  0 


If  we  substitute  this  into  equations  4.2-29  we  find 


0 


CMfl,M  “  RMfl  '  Q 

^-M-1,M  =  ^RMfl  "  Q) '  since  R_M  =  RJ, 


“  eWH,M 


Q .£ .D . 


The  left  part  of  the  flow  diagram  in  Figure  5.1-6 
shows  the  steps  involved  in  the  recursive  computations. 
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4.232  Wiener-Masani  Projections 

The  projection  technique  for  factoring  a  matrix 
valued  autocorrelation  involves  the  theory  of  linear 
algebra.  It  will  be  convenient  for  our  development  to 
consider  vectors  of  matrices  rather  than  polynomials  with 

* 

matrix  coefficients,  i.e.,  we  will  work  in  the  time  domain 
rather  than  in  the  z- transform  domain. 

Let  us  begin  by  defining  the  elements  A.  of 
a  complete  subset  M  of  the  linear  space  of  vectors  of 
matrices  S 

*  V  At+1'  **•'  At+n' 

That  is,  each  element  of  is  a  time  shifted 

reproduction  of  the  minimum  delay  operator  A Q 
(Robinson,  p.  75,  1962). 

We  shall  also  define  an  inner  product 

(Ai-  V  -  As+j-i  a;  (“-2-32) 

■  Ri-J 

This  definition  conforms  to  all  the  requirements  for  a 
linear  product: 

s )  (Ai 9  Aj ) 1  =  (Aj  *  Ai )  * 

b)  (Ai>  Ax)  >  0  if  A^o, 

c)  (a  A^'  Aj)  =  a(A^,  Aj)  where  a  is  a  scalar, 

d)  (Ai  +  Aj  *  Ak)  =  (Ai»  Ak)  +  (Aj  +  A^)  • 
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Spaces  which  are  linear,  complete,  and  contain  an  inner 
product  are  called  closed  linear  manifolds .  Neumann  (p.  51 > 
1950)  states  the  following  definition  and  theorem 
concerning  projections*. 

Definition  4 .2-2 .  If  M  is  a  closed  linear  manifold 
in  S,  if  B  £  S,  and  if  B  =  +  B^,  where  B^  £  M 

and  Bg  €  -M  ,  Then  B-^  is  called  the  projection  of 
B  on  M  and  the  operation  of  projecting  B  on  M  is 
denoted  by  PM  B  =  . 

Theorem  4.2-2.  A  necessary  and  sufficient  condition  that 
an  operator  E  be  a  projection  P  is  that 

a)  E  is  single  valued,  linear  with  domain 
in  S  , 

b)  (E  A^,  Aj)  =  (Ai,  E  Aj)  for  every 
and  Aj  in  S  , 

c)  EE  =s  E  . 

M  is  uniquely  defined  by  E  . 


Finally  Neumann  (1950)  states  the  crucial 
projection  theorem. 

Theorem  4.2-3.  (Projection) .  If  E^  =  P^  and  E^  =  P^, 
then  the  sequence  of  operators  E^,  E,^,  E^E^E^,  E2E1E2E1 
has  a  limit  E  ;  the  sequence  E2,  E-^Eg,  E2E1E2,  E1E2E1E2 


has  the  same  limmlt  E  ;  and  E  =  P 


Ml«2  1 
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Wiener  and  Masanl  (p.  106,  1938)  state  the 
following  corollary. 

Corollary  4.2-3.  If  P  is  the  projection  on  Mi 
then 

P  »  I  "  El  "  Eg  +  E1E2  +  E2E1  -  ... 

The  convergence  being  in  the  strong  sense. 


Wiener  and  Masani  (1958)  then  give  a  lengthy 
development  to  generalize  this  equation  to  include  an 
infinite  number  of  projections.  They  find 


P 


I  -  I  E.  +  E  E.E . 
i=l  1  1, j=l  1  J 


i, J,k-1 


EiE,iEk + 


(4.2-33) 


where  the  projection  operator  is  defined  by  the  inner 
product 


The  normalization  that  they  use  to  insure 
convergence  is  the  requirement  that  i  *  ^  • 

To  make  this  normalization,  we  factor  each  of  the  diagonal 
terms  such  that 

(  R(z)  )ifl  »  ai(z)  rL  a1#(z)  i  -  1,  ...,  n 
where  a^z)  is  normalized  so  that  its  constant  term  is 


equal  to  one. 

Now  let 

a(z) 

=  diag 

[a1(z),  a2(z),  ...,  an(z)J 

y/r 

=  diag 
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then  the  normalized  autocorrelation  is 

R(z)  =  v/r*^a  ^(z)  R(z)  a^fzjv/r’  * 

Our  problem  then,  is  to  find  the  vector,  A  , 
that  is  orthogonal  to  each  of  the  t  =  1,  . 

Thus  we  substitute  A  into  the  projection  sequence  4.2-33: 


FA  =  I  A  -  I  E.A  +  £  2  E.E.F.JA  + 

1=1  1  i,J=l  1  ^  i, j,k=l  1  J  K 


=  I  a  -  £  (a,a^ )A^  +  r  (A,Aj) (Aj,a1)a1 


-  £  (AjA^)  (Aj^Aj)  (A.j*A^)A^ 


eo  f  to  to  n 

A  +  E  -R  .  +  E  R  1  .  -  £  ft  .ft .  .ft  .  +  ...  A. 

1=1  L  1  J=1  J  1  1  J,lo=l  _k  •*  k  1_J  J-1 


Therefore  the  orthogonal  operator 

A°o  ~  ^  1«°  *  *  ••• 

is  given  in  terms  of  the  autocorrelation  only: 


-  1n 


^1-  =  -R-l  +  £  R-J  Vi 


£  ftj,  ft.  b-  ft.  .  + 
k.J-1  k  -1_k 


For  computational  purposes  we  define  a  vector 

1.U  A.  v. 

Cl  j  which  represents  the  J^  projection  of  the  1— 
term  of  A ^  .  Clearly  then,  the  first  projection  is 

*i,0  =  Si 

t .  ,  =  -R  ,  i  *  I#  2,  ... 

1,1  -1 
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and  the  (j  +  1}—  projection  is 

Ci,J+i  =  "C1,J  Ri  "  C2,j  Ri-1  "  "•  "  Ci+1  R0  "  *** 

*i,J+l  =  *i,J  +  Ci,J+l  i  »  1,  2,  ... 

The  iteration  is  continued  until  t,  .  becomes  smaller 

i » J 

than  some  given  value. 

The  only  problem  that  remainB  is  that  of  scaling 

A  nf—  1 

so  that  is  represents  the  inverse  G  (see  equation 
4.2-31).  We  have 

r(z)  *  Kl(z)  A0  ^ii(z) 

and  R(z)  =  £(z)  AQ  £#(z) 

from  the  spectral  factorization  of  a  nc-.-degenerate 
autocorrelation.  But  we  had 

R(z)  *  v/r-*la_1(z)  R(z)  a^fz)  Sr"1 

£(z)  AqA^  ^*(z)  =  a(z  )\/r~  (zJAqA^^1  (z)v'r  a#(z) 

and  ^(z)  »  \/r  Am(z)  v/r*1  a_1(z) 

This  development  is  intended  to  be  a  quick 
summary  of  the  projection  technique.  It  is  by  no  means 
rigorous.  The  step  from  the  Neumann  theorem  to  the  actual 
projection  definitions  that  converge  is  certainly  not 
immediate.  One  must  either  follow  the  path  that  Wiener 
and  Masanl  (1958)  did,  or  generalize  Neumann's  scalar 
theorems  to  matrix  space. 
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4 . 3  Multi -Dimensional  Autocorrelation 


Since  the  properties  of  matrix-valued  multi¬ 
dimensional  autocorrelations  follow  directly  from  those 
of  scalar  autocorrelations  we  will  li*itthe  discussion 
here  to  scalar  values. 

The  properties  of  multi-dimensional  auto¬ 
correlations  have  the  same  features  that  we  observed 
for  one-dimensional  correlations. 

Let  R (z )  =  R(z,  z.,  z^)  be  a 

(k+1) -dimensional  scalar  autocorrelation  function,  then 

a)  R(z)  is  Centro-symmetric;  that  is, 
R(z,z1, . . .,zk)  =  R(l/z,  l/z^ , *  *  * ,  l/zk)  . 

b)  R(z)  is  non-negative  definite  on  the 
unit  hyper-circle  |z  z 1  ...  zk|  =  1  . 

4.31  Mapping  into  One-Dimensional  Representation 

Perhaps  the  most  important  thing  that  we  will 
establish  here  is  the  mapping  of  multi-dimensional 
autocorrelations  into  a  matrix  representation.  We  will 
begin  by  making  the  transformation  in  terms  of  the 
mapping  operators  defined  in  Section  3*^3  and  then 
proceed  to  direct  transformation  from  multi-dimensional 
to  matrix  valued  autocorrelations . 
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The  napping  is  illustrated  in  Figure  4.3-1.  The 
configuration  of  the  autocorrelation  matrix  is  the  same 
that  we  would  obtain  if  we  had  treated  the  process,  x(z), 
as  three  separate  wavelets,  xi(z),  Xgfz),  and  x^,(s) 
and  defined  the  autocorrelations  as 


Xj(z) 

x1#(z),  x2#(z),  x3#(z)l  = 

»  - 
r0(z)  rx(z)  r2(z) 

^(z) 

J 

r_l(z)  rQ(z)  r3(z) 

x3(z) 

r_2 (z)  *  x(z)  rQ(z) 

•  m 

However,  for  each  spatial  lag  we  must  take  the  sum  of  the 
correlations  of  all  the  wavelets  that  overlap  at  that  lag. 
The  matrices  shown  in  Figure  4.3-1  fulfill  thiB  require¬ 
ment,  They  are  defined  formally  by 

MAP2 (x,  x (z ) )  *  MAP? (x,  x(z))»  *  R(z)  . 

Frequently,  in  practical  applications  we  are 
presented  with  the  multi-dimensional  autocorrelation  r(z) 
and  we  wish  to  map  it  directly  into  the  matrix  represen¬ 
tation  K(z)  .  The  procedure  here  is  very  similar  to  that 
taken  above.  We  map  the  spatial  positions  of  the  process 
into  a  vector  and  form  the  symbolic  product 


V 

[*1  X2  *N] 

3 

ri.a 

**•  rl,N 

x2 

• 

• 

r.,i 

• 

• 

r2,2 

*•*  r2,N 

•  • 

•  * 

• 

-XN  - 

• 

_rH,l 

r 

N,c 

•  • 

**•  rN,N. 

1  J  * 

-L‘“  J 


T 


‘t  1  * 


© 


t #  i^» ig 


=  autocorrelation 


L 


z-Transform  in  preferred  direction 
and  order  the  indices. 

x1(z)  x3(z)  x1#(z) 

x3*(z)  r 1 

x2(z) 

x2#(z) 

r-2<z>  r0 

r_3(z)  r- 

i 

Map  wavelets 

into  matrix  notation. 

r* 

.  n  r“ 

1 

"t  /• 

Mz;  r3(z) 


Xj  Xj  0  x3  0  <T 

X 

M 

* 

o 

o 

_ 1 

~r0  ll  r2~ 

0  X-.  x„  0  x-j  0 

1  ^  3 

! 

V  xi»  0 

2 

r-l  ro  r3 

[0  0  0  x!  *2  x3j 

0  Xj,  0 

*3«  0  xi» 

®  x3*  x2» 

r-?  r-3  r0 

0  0  x3. 

Figure  4.3  -  1:  Mapping  of  a  multi -dimensional 
autocorrelation  into  a  matrix  representation. 


124 


define  the 


Then  the  subscripts  of  each  element  r 
spatial  separation  for  that  term.  If  the  spatial  process 
is  stationary.,  that  la,  if  r^  ^  depends  only  upon  the 
spatial  separation  between  position  i  and  position  J, 
then  R(z)  Is  symmetric.  This  leads  to  a  number  of 
simpllf 1 cations . 

4.32  Methods  of  Factorization 

Except  for  one  method  of  factorization,  all  of 
the  techniques  that  are  used  are  made  in  terms  of  the 
matrix  mapping  of  the  autocorrelation.  The  fact  that  the 
autocorrelation  matrix  is  symmetrical  may  lead  to  some 
important  simplifications  in  some  cases.  For  example,  in 
the  elementary  autocorrelation  matrix  factorization, 
synsnetricallity  forces  the  algorithm  to  give  a  one-sided  t 
wavelet.  Also,  in  the  least- squares  approximate  technique 
the  operator  >fM(z)  -  $M(z). 

It  is  instructive  to  consider  the  meaning  of 
the  minimum  delay  wavelets  that  one  obtains  from  the 
matrix  factorization.  Each  row  of  the  matrix  A(z)  will 
be  a  vector  representation  of  a  spatial  minimum  phase 
wavelet.  This  vector  representation  is  the  same  as  that 
used  for  mapping  the  original  process.  The  origin  of 
the  wavelet  is  located  at  the  spatial  position  corresponding 
to  the  diagonal  term  in  the  matrix.  Thus  the  autocorrelation 
of  a  spatial  process  having  n  lattice  points  will 
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produce  n  minimum  phase  wavelets;  and  eacn  wavelet  well 
have  its  origin  at  a  different  lattice  point. 

Occasionally  in  physical  problems  we  know  that 
the  factorization  should  have  zero-phase,  i.e.,  snould  be 
symmetrical  in  all  directions.  For  this  case  we  may  proceed 
as  in  the  one -dimensional  scalar  case.  Thus  we  need  only 
evaluate  the  expression 

a  (e*1-)  =  . 

That  is,  we  find  the  cosine-transform  of  r^  ,  take  its 
square  root,  and  retransform  back  to  space-time. 


.1 


5.  LEAST-SQUARE  FILTERING  IN  THE  PRESENCE  OP  NOISE 

In  Chapter  4.  a  number  of  techniques  were 
discussed  for  finding  a  minimum  delay  wavelet  from  a  given 
autocorrelation.  Of  the  techniques  discussed,  the  least- 
squares  approximation  was  found  to  be  the  best  method  In 
the  sense  of  computational  efficiency.  In  thi3  chapter, 
the  least-squares  decompositional  method  will  be  extended 
to  include  signal  shaping  (in  addition  to  straight  predic¬ 
tion)  in  the  presence  of  random  noise  with  a  given 
coherency.  This  approach  will  give  an  optimum  linear 
operator  for  a  given  length  and  output  lag. 

The  normal  equations  for  the  one-dlmenslonal 
matrix-valued  process  only  will  be  developed  here.  As  was 
illustrated  In  the  last  two  chapters,  all  other  dimension¬ 
alities  are  but  a  special  case  for  this  representation. 

5.1  Derivation  and  Recursive  Solution  of  the  Normal 

Equations  ‘ 

The  solution  of  the  problem  of  determining  the 
optimum  least-squares  linear  operator  is  based  upon  the 
following  assumptions: 

a)  The  known  n  x  m  matrix-valued  signal  St 
is  the  additive  combination  of  K  uncorrelated  stationary 
random  processes  3,  *  . 

1 ,  w 
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b)  The  n  x  m  matrix-valued  noise  ]if  is  a 
random  process  with  z 
covariance,  E  NJ 

c)  The  observed  random  process,  X^.  ,  is  the 
additive  combination  of  the  signal  and  the  noise 

Xt  "  St  +  Nt  . 

d)  The  observed  random  process  is  convolved 
with  an  undetermined  l  x  n  matrix-valued  wavelet  F1 
i  »  1,  . M  to  obtain  dhe  l  x  m  matrix-valued 
actual  output  Yf  . 

e)  The  l  x  m  matrix-valued  desired  output, 

,  is  the  additive  combination  of  K  independent  desired 

outputs  ft,  .  where  15,  .  is  unccrrelated  with  % ,  .  , 

1,U  1,V  J  ,  V 

1/  J  ,  l.e. 

E  {^i,t+T*j,o)  =  0  • 

5.11  Normal  Equations 

The  linear  least-squares  operator  wavelet  is 
determined  by  requiring  that  the  norm  of  the  difference 
between  the  actual  output  and  the  desired  output  is 
minimum  for  all  time.  That  Is,  we  require  that  be 

mini  .lized,  where 

=  2  ( (  e^.  ) ) 

-  E  ( tr  et  -t)  - 


ere  mean. 


)  • 


<V 


=  0  ,  and  known 

n ,  m 
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et  “  Dt  "  Yt  * 

and  E  stands  for  expected  value. 


To  find  the  minimum,  we  take  the  derivative  of 
with  respect  to  the  coefficients  of  the  wavelet  Fj 
J  »  1,  . ..,  M  and  set  it  equal  to  zero.  Thus 


J  -  1, 


•  •  •  $ 


M 


implies  that  the  error  e^  is  normal  to  the  input  X^j^: 

E  (el  Xl-J+l)  ■  0  4  -  1,  ...»  M  . 

This  orthogonality  was  the  basis  for  the  development  of  the 
Wiener-Masani  projections  (see  Section  4.232).  This  is  alBO 
the  origin  of  the  name  "normal  equations." 


Now,  let  us  expand  the  normal  equations: 

E  {ei  xi-,j+i)  "  0  J  -  i)  •••,  m 

-  e  (  (di  -  pk  xi-k+i>xi-J+i ) 

Also,  we  have  +  N^  and  B<Ni)  •  Thus  the 

normal  equations  have  the  form 

E  {Di  ^Si-J+1  +  Nl-j+l^ 

"  ^  pk  ^Si-k+l  +  Ni-te*-l^Si-j+l  +  Ni-j+l^)-  0 
E(^l  Pk  (Sl-k+l  Si-J+1  +  Ni-k+l  Ni-j+l^“  E  (Di  Si-j+]) 
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E  (SJ-k  so)'L  E  (NJ-k  No)  ■  E(DJ-lSi) 

1,  M  .  From  the  assuptions  we  see  that 

E  (SJ-k  S0>  ■  £  E<*l.J-kSi.o> 

E(Dj-kso)=  E(T,i,j-kSi,o)  • 

Tnerefore  the  autocorrelation  of  Sfc  is  the  sum  of  the 
autocorrelations  of  3,  .  .  If  we  define  an  auto- 

1  j  J 

correlation 

R1  -  E  (Si  so)  +  E  (Ni  «o) 
and  a  cross-correlation 


then  the  normal  equations  may  be  written  in  the  simple 
form 

m  Fi  RJ-k  “  “j-1  *  J  =  1 . M  •  (5a"1) 


M 

F  Fk 
te>l 

for  J  = 


We  may  also  obtain  a  simple  form  for  the 
expected  error  : 


tr 

tr 

tr 


M 

“  1  Fl 

k*l  1 


since 

Xi-k+l 


Xi-k+l^  ^ 

e^  is  noririal  to 
for  k  =  1, . . M 
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'•o  ■  tr  E  (Di  D1  '  ii  xi-k+i  Di) 

“1  •  1  '«  °w>  .  B.M) 

The  normal  equations  that  were  obtained  above  are 
very  closely  related  to  the  simultaneous  equations  defined 
in  Section  4.231  for  A.  „ 

1|N 

M 

k*0  j  *»  0,  *,,,  K 

where  AQ*M  =  I  ,  and  Cj  M  *  0  ,  J  -  1,  ...,  M  .  For  if 
there  is  no  noise,  =  0  for  all  i  ,  and  if  we  make  the 

desired  output  equal  to  the  input  one  lag  ahead  in  time 
(that  is,  we  ask  the  filter  to  predict  the  next  value  of 
the  process)  then  the  n  x  n  filter  wavelet  i  -  1,..,M 

is  identical  to  the  wavelet  -A^  M  i  »  1,  ...,  M  and  the 
equations  above  for  J  *  1,  ...,  M  are  the  same  as 
equation  5.1-1.  Also,  the  equation  for  the  expected  error 
m  is  the  same  as  equation  5.1-2. 

The  filter  -A^  ^  i  *  1,  ...,  M  is  called  the 
least-squares  approximate  prediction  filter  with  unit 
prediction  distance.  The  difference  between  the  actual 
value  and  the  predicted  value  is  the  output  of  Aa  u 
1  =  0,  ...,  M  .  Thus  it  is  called  the  prediction  error 
filter  or  the  foresight  error  filter. 
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On  the  other  hand,  had  we  first  reversed  the 
direction  of  time  and  then  solved  for  the  prediction  error 
filter  as  above,  we  world  have  obtained  the  filter 

5i,M  1  *  °*  M  : 


M 

ufo  Rk-J  ■  ?J,M 


J=0,  * . . ,  M 


where  BQ  M  *  I  and  Ej  M  =  0  J  *  1,  ...»  M  .  Since  this 
filter  actually  predicts  the  past  from  future  values  of  the 
prccess,  it  is  called  the  hindsight  error  filter. 


Example  5.1-1*  Let  us  consider  a  symbolic  two-dimensional 
scalar  problem  of  signal  shaping  in  the  presence  of  noise. 
The  problem  is  specified  by  the  arrays  illustrated  in 
Figure  5.1-1.  That  is,  we  wish  to  design  a  filter  which 
produces  the  desired  output  array  when  convolved  with  the 
signal  array  and  which  produces  a  zero  output  array  when 
convolved  with  the  noise  array. 


The  actual  design  of  the  filter  is  based  on  the 
two-dimensional  autocorrelations  of  the  signal  and  of  the 
noise,  and  on  the  two-dimensional  crosscorrelation  of  the 
desired  output  with  the  signal.  These  arrays  are  shown 
in  Figures  5.1-2  and  -3. 


These  arrays  are  then  mapped  into  matrix 
notation  (see  Cnapter  4.31)  and  substituted  into  the 
normal  equations.  The  recursive  procedure  described  in 
the  next  section  was  then  followed  to  compute  the  20  x  20 
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Autocorrelation  of  the  Noise  Array 
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Figure  5*1  “  2:  Complete  2-dimensional  autocorrelations 
of  the  signal  array  and  of  the  noise  array. 
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term  filter  shown  in  Figure  5*1-4.  This  recursion  was 
carried  out  by  first  extending  the  length  of  the  filter  in 
a  positive  direction  and  then  extending  it  in  a  negative 
direction  alternately  until  the  final  length  was  reached. 
After  each  extension,  the  expected  error  was  computed. 

A  plot  of  these  values  is  also  shown  in  Figure  5*1-4 • 

Finally,  to  gain  a  visual  idea  of  the  quality 
of  the  filter,  it  was  convolved  (in  two  dimensions)  with 
*he  signal  array  and  with  the  noise  array  to  produce  the 
two  output  arrays  shown  in  Figure  5*1-5*  Only  the  center 
portions  of  the  convolutions  are  shown. 
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Figure  5.1  -  3:  Complete  crosscorrelation  of  the  signal 
array  with  the  desired  output  array. 
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Figure  5.1  -  **:  Optimum  20-row  x  2C-coiumn  filter  com¬ 
puted  from  the  correlations  in  Figures  5.1“2  and 
The  expected  error  is  for  20  x  1#  20  x  2, 
...»  20  x  20  term  filters. 
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Convolution  of  Filter  with  Sign*!  Array 
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Convolution  of  Filter  with  Noise  Array 
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Figure  5.1  -  5s  Actual  outputs  of  filter  when  it  is  con¬ 
volved  with  the  signal  array  and  with  the  noise  array. 
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5.12  ,  5.  121 


5.12  Recursive  Computation  Algorithm 

The  recursive  scheme  that  was  presented  In 
section  ^.231  can  be  expanded  to  apply  to  the  extension 
of  the  length  of  the  wavelet  .  In  addition,  Simpson 
(Simpson,  et  al,,  1963)  has  proposed  a  similar  recursive 
scheme  to  shift  the  lag  between  the  desired  output 
and  the  input  X<  .  Thus  this  recursion  allows  and 
efficient  search  for  the  optimum  lag. 


5 . 12 1  Extension  of  filter  length 


The  recursive  algoritnm  will  be  stated  in  terms 
of  the  z-transform.  The  normal  equations  for  FM(z)  roay 
be  written  as 


Fm(z)  R(z)  =  eH(z) 


where  ^  =  ^J*l  ^  ^ *  •••»  M  • 


If  we  weight  B^l/z)  and  add  it  to  FM(z ) 


we  find 


[pH(z)  +  z^1  Kp^^d/z)]  R(z)  -  eM(z)  ♦  z*1  Kp^K^z) 
Since  =0  i  =  -1,  . . . ,  -M  ,  if  we  choose  Kj^m 


such  that 


*Mfl,M  +  KP,M^0,M  *  GM  ' 


then  we  obtain  the  recursive  relationship 
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Figure  5.1  -  6:  Recursion  to  extend  the  length  of  an 
optimum  least-squares  filter.  The  numbers  on  the 
boxes  illustrate  a  possible  computational  sequence. 
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Vl<2)  ■  PM<Z>  +  zMfl  S.H  ' 

This  recursion  is  Illustrated  in  detail  in  Figure  5.1-6. 

Similarly,  if  we  had  wished  to  extend  F^(z)  in 
the  other  direction,  we  would  weight  A^{z)  and  add  it  to 

PM(z) 

PM<Z>  +  VZ>  *  eM<Z>  +  «?,H  VZ>  ' 

A 

where  we  choose  K„  M  such  that 

F,  M 

£0,M  +  *F,M  e0,M  =  G-1 

and  we  have  extended  the  length  of  ( z )  in  the  negative 
direction  according  to 

F«fl<z>  “  z  [PM(z)  +  *F,M  VZ)1  ’ 

5.122  Shift  of  output  origin 

This  recursion  is  slightly  more  complicated  than 
the  simple  length  extension.  Our  objective  is  to  find  a 
filter  such  that 

F^fz)  R  (?)  =  et1J(z) 

where  =■  Gj_2  j  ■  1,  ....  M  given  FM(z)  as 

described  in  the  last  section.  The  first  step  is  to  shift 

fm(z)  and  B  (1/z)  right  one  unit  and  subtract 

M  M-l 

zK<z)  “  zM  FK,M  ^M-l^1^2^]  R^Z-  =  *VZ^ 
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where  ^  M  =  Gi_2  i  *  2,  . ..,  M  .  Thus  the  error  that 
is  Introduced  by  shifting  PM(z )  Is  compensated  for, 

1  =  2,  M,  by  subtracting  multiplied  on  the 

left  by  ?M  M  .  Now  we  add  a  weighted  version  of 

to  adjust  the  value  of  Jf,  M  to  be  G  ,  : 

i,n  “1 

2  [PM^Z ^  “  Z  PM,M  5M-1^1//z^  *  *M-1^]  * 

(Z )  +  ^FJM  ^M-l^ 

Thus  If  we  i-  it 

*1,M  +  C0,M-1  "  °-l 

we  find  that 

PM  )(z)  -  z  [PM(z)  -  zM  pmsM  +  ^JM^M-l^J  * 

Similarly,  we  could  have  chosen  to  left  shift  the 
cross-correlation  to  solve  for  F^”^(z) 

t 

H(z)  -  ejt"1^) 

*  (-1) 

where  J  ■  1,  ...»  M.  This  filter  is  given  by 

-  1/z  [p„(z)  -  V1  #M  z  4^(2)  + 

V.n  Vi*1/*)] 

where  yj  -  Vl,K-l  ”  °M  *  *nd 
*  ^P2,M  “  P1,M  41,M-1^  R-1  * 

* 

+  (fh,m  "  vx,n  ^w-i,ii-i^  R-ir 

*t 
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5.2  Computational  Properties 

The  normal  equations  derived  in  the  first  part 
of  this  section  will  certainly  provide  the  optimum  linear 
least-squares  filter  for  a  given  specification  of  a  problem. 
The  success  of  the  application  of  such  filters  depends  to 
a  large  extent  upon  the  design  of  the  specifications  so 
that  a  small  expected  error  is  obtained.  The  rest  of  this 
chapter  will  be  devoted  to  examining  some  qualitative 
design  criteria  for  scalar  least-squares  problems.  Most 
of  the  conclusions  may  be  easily  extended  to  multi-variable 
problems. 


In  most  least-squares  filter  problems  there  are 
at  least  three  parameters  that  are  left  to  the  discretion 
of  the  designer: 

1)  The  inclusion  of  an  arbitrary  amount  of 
white  noise  to  represent  uncertainty  in  the 
design  criteria. 

2)  The  specification  of  the  shape  of  the 
desired  output  process. 

3)  The  specification  of  the  lag  between  the 
desired  output  and  the  input  process. 

A  set  of  computational  experiments  were  performed 
to  test  the  effects  of  the  variation  of  these  and  other 
parameters.  In  each  experiment,  after  an  autocorrelation 
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series  and  a  crosscorrelation  series  were  specified,  the 
expected  errors  for  all  filters  with  lengths  less  than 
some  maximum  (50  points  in  these  examples),  and  with  all 
relevent  lags,  were  computed.  The  array  of  expected 
errors  were  normalized  between  1.  and  0.  and  were 
contoured  in  terms  of  decibels  (10  los^o  ®  where  ?.  is 
the  expected  error) .  Thus  each  experiment  provides  a 
complete  test  of  the  effect  of  the  output  lag. 

Nearly  all  of  the  experiments  were  performed 
using  the  wavelets  shown  in  Figure  5*2-1.  The  mixed 
delay  and  minimum  delay  wavelets  showu  have  the  same 
amplitude  spectrum.  The  mixed  delay  wavelet  was  chosen 
with  a  large  dynamic  range  so  that  the  expected  error 
plots  would  have  definite  character. 


Amplitude  Spectrum 


0  T 


K*  ;ure  -  Is  'line  domain,  z-transform,  and  amplitude 
spectrum  of  a  mixed  delay  wavelet  and  its  equivalent 
minimum  delay  form. 


3.2 1 


3.21  White  Noise 

The  first  experiment  involved  signal  shaping 
in  the  presence  of  white  noise.  Since  the  autocorrelation 
of  white  noise  is  a  spike  at  zero  lag  and  zeros  elsewhere, 
white  noise  is  included  by  adding  a  constant  to  the  center 
term  of  the  autocorrelation. 

The  filters  for  which  the  expected  errors  were 
computed  were  asked  to  compress  the  mixed  delay  wavelet 
shown  in  Figure  5.2-1  into  a  spike  while  rejecting  varying 
amounts  of  white  noise.  The  resulting  expected  error 
arrays  are  shown  in  Figure  5*2-2.  The  primary  point  to 
notice  here  is  that  as  the  relative  power  of  the  white 
noise  is  increased,  the  filter  shape  becomes  stable  for 
shorter  lengths  but  the  expected  error  approaches  a 
constant  greater  than  zero.  This  constant  is  related  to 
the  relative  powers  of  the  white  noise  and  the  wavelet. 

The  position  of  the  minimum  relative  to  the  output  lag  is 
also  somewhat  dependent  upon  the  amount  of  noise  added. 
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5.22 


5.22  Delay  Properties  of  the  Input  Wavelet 

Although  it  is  usually  not  readily  alterable*  the 
delay  of  the  input  wavelet  will  certainly  affect  the  design 
of  the  filter  problem.  This  experiment  illustrates  the 
effect.  For  each  case  the  least-squares  filter  was  asked 
to  compress  the  energy  of  the  input  wavelet  into  a  spike. 
The  contours  of  the  expected  error  arrays  are  shown  in 
Figure  5.2~3.  The  first  plot  is  for  the  minimum  delay 
wavelet  shown  in  Figure  5.2-1.  The  second  plot  is  for  the 
mixed  delay  wavelet,  and  the  third  plot  is  for  a  maximum 
delay  wavelet  obtained  by  taking  the  time  reverse  of  the 
minimum  delay  wavelet. 

Clearly  the  location  of  the  minimum  is  strongly 
affected  by  the  delay  of  the  input  wavelet.  The  relation¬ 
ship  between  the  position  of  the  expected  error  minimum 
and  the  delay  is  not  simple.  A  longer  set  of  experiments 
showed  that  this  position  of  the  minimum  is  most  highly 
depen  lent  upon  the  zeros  that  are  far  from  the  unit  circle 
and  least  dependent  upon  seros  close  to  the  v  .it  circle. 
That  is,  the  change  in  position  of  the  minimum  when  a 
distant  root  is  reflected  about  the  unit  circle  is  much 
greater  than  when  a  root  with  nearly  unit  magnitude  is 
reflected.  It  is  possible  to  say,  though,  that  for  all 
the  possible  minimum  positions  for  wavelets  with  a  given 
amplitude  spectrum,  vhe  minimum  delay  "‘velet  will  have 
the  position  closest  to  the  left  side  of  the  plot,  and 
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Figure  b.2  -  3:  Contours  (in  decibels)  of  the  expected 
error  arrays  of  spiking  filters  for  3  wavelets  that 
have  the  same  amplitude  spectrum. 


5.23 


the  maximum  delay  wavelet  will  have  the  position  closest  to 
the  right  side.  As  the  expected  error  at  these  optimum 
lags  approach  zero  the  optimum  position  will  approach  lag 
one  for  a  minimum  delay  wavelet. 


5.23  Desired  Output  Spectrum 

The  last  set  of  experiments  tests  the  effect 
on  the  expected  error  of  the  shape  of  the  desired  output 
wavelet  relative  to  that  of  the  input  wavelet.  The 
definition  of  such  an  experiment  is  necessarily  very 
vague  since  it  is  difficult  to  define  a  measure  of  relative 
shape.  For  the  examples  shown  in  Figure  5.2-4  the 
amplitude  spectrum  of  the  desired  output  was  varied  relative 
to  that  of  the  input  wavelet. 

The  input  wavelet  for  all  of  the  cases  was  the 
mixed  delay  wavelet  shown  in  Figure  5.2-1.  The  desired 
output  for  the  first  case  had  the  same  amplitude  spectrum 
as  the  input  wavelet  but  was  maximum  delay.  The  desired 
output  for  the  second  case  was  a  spike,  that  is,  it  had 
a  unit  amplitude  spectrum.  For  the  last  case  the  desired 
output  was  a  20-term  minimum  delay  wavelet  that  had  an 
amplitude  spectrum  that  was  approximately  reciprocal  to 
the  amplitude  spectrum  of  the  input  case. 

This  set  of  examples  verifies  the  conclusion 
that  one  would  make  intuitively,  i.e.  that  the  closer 
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Figure  -  4:  Contours  (in  uecibels)  of  the  expected 

error  arrays  of  optimum  shapln  filters  for  whlcn  the 
desired  output  nau  the  same#  a  flat#  and  ai-  Inverse 
amplitude  spectrum  relative  to  the  signal  wavelet. 


the  spectrum  of  the  desired  output  wavelet  is  to  that  of 
the  input  wavelet,  the  better  the  filter  is  able  to  perform. 

5.24  Output  Lag  and  Filter  Length 

All  of  the  examples  illustrated  in  Figures  5.2-2, 
-3*  and  -4  illustrate  properties  of  expected  error  verses 
output  lag  and  filter  length.  We  can  draw  the  following 
conclusion:  The  expected  error  is  a  non-increasing  function 
of  filter  length  for  any  particular  lag  of  the  desired  output 
relative  to  the  input.  The  value  of  the  expected  error 
levels  out  to  some  value  which  depends  on  the  output  lag. 
Indeed  the  expected  error  is  strongly  dependent  upon  the 
output  lag  and  may  vary  quite  rapidly  for  small  changes  in 
the  output  lag.  Generally,  as  the  length  of  the  filter 
becomes  long  with  respect  to  the  lengths  of  the  input  and 
desired  output  wavelets,  the  expected  error  curve,  plotted 
with  respect  to  the  output  lag,  has  one  relative  minimum. 

The  discussion  above  may  be  sharpened  somewhat 
by  examining  the  results  of  Claerbout  and  Robinson  (1964). 
They  showed  that  the  sum  of  the  expected  errors  for  al1 
possible  lags  of  the  desired  output  relative  to  the  input 
Is  Independent  of  the  length  of  the  filter  if  the  problem 
Involves  no  noise  suppression.  Thus  as  the  filter  becomes 
longer  the  total  expected  error  will  be  spread  over  a 
greater  region.  Consequently,  there  must  be  some  lag  for 
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vhic.v  th i  expected  error  approaches  zero  at  least  as  1/to 
where  in  -s  the  length  of  the  filter. 

The  examples  illustrated  in  Figures  5.2-2,  -3, 
and  -4  show  that  the  performance  of  an  optimum  filter  is 
strongly  dependent  upon  the  design  of  the  input  parameters. 
Thus  the  successful  application  of  least-squares  techniques 
should  include  all  possible  means  of  optimizing  the  design 
of  tne  problem  before  finding  the  optimum  filter  for  a 
particular  case. 


O.  ,  6.1 

6.  GEOFHYSICmL  applications 

The  mathematics  developed  in  the  preceeding  three 
chapters  is  general  in  its  applicability  to  geophysical,  as 
well  as  non-geophysical,  problems.  Since  the  range  of 
possible  applications  is  great,  this  chapter  will  be  restric¬ 
ted  to  two  examples  for  which  the  setup  of  the  problems  do 
nov  require  long  derivations. 

The  first  example  is  an  illustration  of  the  use¬ 
fulness  of  the  zero-phase  multi-dimensional  factorization. 

The  second  example  is  of  a  two-dimensional  least-squares 
problem  that  is  derived  with  a  slightly  different  initial 
criteria  than  that  used  in  Chapter  5.  The  end  result, 
however,  is  a  set  of  normal  equations  for  which  the  recur¬ 
sive  procedures  given  in  Chapters  4  and  5  are  applicable. 

6 . 1  World-Wide  Average  Gravity  Anomalies 

The  Army  Map  Service  (1959)  has  reported  on  a 
statistical  analysis  of  available  world-wide  gravity  data. 
Part  of  the  results  reported  were  numerical  estimates  of 
the  average  gravity  covariance  for  continental  regions  and 
for  oceanic  regions.  These  estimates  are  shown  in  Figure 
6.1-1.  The  curvas  are  an  eye-ball  smoothing  of  the  data. 

If  we  follow  the  assuption  of  the  Anny  Map  Service 
that  the  shape  of  gravity  anomalies  are  stationary  with 
respect  to  azimuth,  then  the  curves  illustrated  would 
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represent  one-half  of  the  cross-section  of  the  two- 
dimensional  autocorrelation  of  an  average  gravity  anomaly. 
Clearly  we  should  be  able  to  apply  the  zero-phase  (see 
Section  4.32)  factorization  technique  to  determine  the 
average  synmetrlcal  gravity  anomaly.  This  was  carried  out 
as  follows: 

1.  A  two-dimensional  autocorrelation  array  was 
made  by  sampling  the  curves  in  Figure  6.1-1  at  1°  intervals. 

2.  The  two-dimensional  cosine  transform  was 

computed. 

3.  Since  the  sample  arrays  are  covarience 
estimates,  they  may  not  represent  possible  autocorrelation 
functions,  that  is,  some  values  of  the  cosine  transform 
may  be  negative.  In  this  example,  a  few  small  values  for 
the  oceanic  samples  were  negative.  These  values  were 
arbitrarily  set  positive.  (As  a  check,  these  values  were 
also  set  to  zero.  The  resultant  solution  was  modified 
only  slightly.) 

4.  The  square  root  of  each  value  in  the  trans¬ 
form  was  taken. 

b*  These  values  were  than  inverse  transformed 
to  obtain  an  average  gravity  anomaly. 

Cross-sections  of  the  resultant  average  gravity 
arrays  are  plotted  in  Figure  6.1-2.  A  preliminary  inter¬ 
pretation  was  made  by  determining  the  depth  of  a  point 
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or  continental  regions 


source  that  would  give  the  same  gravity  values  at  0*  ana 
at  1°.  Such  sources  would  have  a  depth  of  -^60  km  for  the 
continental  regions  and  ~80  km  for  the  oceanic  regions. 

The  anomalous  masses  would  be  ~10£  larger  in  the  oceanic 
regions  than  in  the  continental  regions. 

The  point  source  interpretation  does  not  account 
for  the  large  values  of  the  oceanic  anomaly  out  to  7°. 

The  total  mass  of  an  anomalous  body  is  proportional  to  the 
volumn  under  the  two-dimensional  gravity  anomaly.  Using 
this  property  as  a  criterion  for  comparing  masses ,  we  see 
that  the  oceanic  anomalous  masses  are  approximately  3.5 
times  as  large  as  the  continental  masses. 
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The  idea  of  band-pass  or  band-reject  filtering 
is  perhaps  one  of  the  oldest  concepts  in  one-dimensional, 
scalar  filter  design.  This  concept  can  be  readily  ex¬ 
tended  t?  two-dimensional  filtering  by  defining  band¬ 
pass  or  band-reject  areas  in  the  two-dimensional  Pourier- 
transform  domain  (sometimes  called  the  m-k  plane).  One 
very  simple  form  of  such  filters  has  wide  application  in 
exploration  seismology  for  descrimination  between  plane 
waves  on  the  basis  of  the  direction  of  arrival  of  the  wave 
(see  Fall  and  Grau,  1963;  or  Embree,  et  al,,  1963). 

In  this  section  we  will  consider  the  general 
problem  of  optimally  designing  a  discrete  two-dimensional 
filter  from  band-pass  and  band-reject  area  specification 
in  the  m-k  plane. 

Let  us  first  consider  the  two-dimensional  band¬ 
pass  filter.  The  criterion  for  designing  this  filter  is 
that  we  wish  to  minimize  the  square  of  the  difference 
between  the  desired  band-pass  configuration  and  the  actual 
Fourler-transform  of  the  filter.  Thus,  if  we  let  f j 
J  *  1,  ,..,  n  represent  a  discrete  filter  with  lattice 
points  at  (Xj,  t^)  J  -  1,  ...,  n,  and  d(u>,t)  represent 
the  desired  band-pass  conf iguration,  then  the  error  to  be 
minimized  is 


6.2 


ft 


J  f  f  f  fj  eit,"tj‘kx.1)  -  d(m,k)J  dm  dk  (6.2-1) 


where  a c  =  a a  .  Taking  the  derivative  of  ft  with  respect 
to  f y  we  find  a  minimum  given  by 


2ft 

“FT 


0 


J 


J  =  1|  •  «  .  i  n 


or 


0  = 


r  r  ex<»w  • 

^  -T7 


fj  -  d(m.k)l  d~  dk 


J  a  li  *••»  n,  f,  —  ij  •  •  •  j  n  (6.2-2) 
which  leads  to  the  particularly  simple  result: 

fj  =  -^[  f  e1  («,,tj“lcCj)  dKk)  du»  dk 

J  4TT  '-TT  '-TT 

Thus  the  filter  coefficients  are  determined  by  the 
Fourier- transform  of  the  desired  band-pass  configuration. 
The  expected  error  ft,  is  given  by  (beginning  from  equation 
6.2-1) 

ft  =  J"1  /n  dK*)  £  Jj  fj  51(<8VkXJ)  -  HTS7PT]  cimdk 

which  gives  (according  to  equation  6.2-2) 


f  f  d2(u),k)  dm  dk 

' -TT  J-TT 


r  f4 
J  ' 
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This  is  reminiscent  of  the  expression  for  the  expected 
error  that  was  found  in  chapter  5.11. 


The  simplicity  of  this  method  of  filter  compu¬ 
tation  makes  it  feasible  to  examine  the  properties  of  such 
filters  in  several  ways.  Figures  6.2-1  and  -2  illustrate 
several  of  the  possible  computations  that  are  useful  in 
studying  the  properties  for  any  given  u>-k  configuration. 
Each  of  the  figures  contains  4  contour  plots  that  represent 

1.  Hie  coefficients  for  a  square  13  x  13  term 

filter. 

2.  The  actual  Fourier  transform  of  the  13  x  13 
term  filter  contoured  in  decibels.  This  plot  is  super¬ 
imposed  upon  a  diagram  of  the  desired  pass  region. 

3.  The  array  defined  by 


e 


J 


d2(t»,k) 


duu  dk 


contoured  at  values  of  e^  =  0.00001,  0.0001,  0.001,  0.01, 
and  0.1  .  The  value  at  any  point  in  thin  array  is  just 
the  value  by  which  the  normalized  expected  error  is  decreased 
by  the  addition  of  that  point  to  the  filter.  Thus  these 
contours  represent  the  optimum  filter  shapes  for  minimizing 
both  the  expected  error  and  the  number  of  filter  coefficients. 
It  is  interesting  to  notice  that  the  maxima  in  the  e^  array 
are  parallel  to  the  boundaries  of  the  pass  region  in  the 
m-k  definition  of  the  filter. 
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Number  of  Traces 


6.2 


•* 


.  Length  of  Traces 

Figure  6.2  -  1;  The  filter  coefficients  for  an  optimum 

least-squares  band-pass  filter;  the  contours  (in  deci¬ 
bels)  of  the  actual  #!-k  plane  reporise  of  the  filter: 

fjpect??  orror s  for  a11  rectangular  fil- 
ters  smaller  than  D 1  x  bl  terms;  and  the  optimum 
filter  shapes  *  or  this  pass-band  configuration. 


i 


Number  of  Traces 


Length  of  Traces 

Figure  6.2  -  2:  The  filter  coefficients  for  an  optimum 

least-squares  band-pass  filter;  the  contours  (in  deci¬ 
bels)  of  the  actual  u>-k  plane  response  of  the  filter; 
the  a  "ray  of  expected  errors  for  all  rectangular  fil¬ 
ters  smaller  than  51  x  51  terms;  and  the  optimum 
filter  shapes  for  this  pass-band  configuration. 


i 


o.2 


4.  The  array  o  normalized  expected  errors  for 
all  possible  rectangular  dimensions  of  the  filters  smaller 
than  51  x  51  terms.  The  center  dotted  line  represents 
the  optimum  length  for  any  given  number  of  traces  while  the 
area  between  the  outer  dotted  lines  include  all  dimensions 
which  are  quite  close  to  the  optimum  rectangular  shape. 

In  nearly  all  physical  problems  the  contribu¬ 
tions  to  the  error  in  some  parts  of  the  m-k  plane  are 
more  important  than  the  contribution  in  other  parts  of  the 
plane.  Let  us  define  a  weighting  function  W(w,k)  that 
expresses  this  importance.  The  larger  W(u>,k)  is  in  any 
given  area,  the  more  important  the  error  is  in  that  area. 
Let  us  Introduce  W(u>,k)  into  the  least-squares  band-pass 
problem: 

ft  *  /  W{tt»,k)  ^  f j  e^V^j)  -  d(»,k)j  2  dm  dk 

Setting  the  derivatives  of  ft  to  zero  as  before,  we  find 

—  =  o 

-  J  f  W(iu,k)  eMu>tt-kxt)  ' 

-rr  tt 

U  fj  ,-M-ytaj)  .  d(„,ic)J  *,  dk 

we  find 
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b.2 


/tt  ,  tt 

J  W(u),!c)  exp  -  k(x^~Xj)J  du>dk 


/  /  «<  tu,k)  d(tu,k)  exp  i^utt^-kxj  diwdk 

-TT  ''-TT 


^  a  I#  •••>  n 


Furthermore,  the  expected  error  becomes 


l  =  J  J  W>,k)  d(u>,k)^£  fj  e“1^u,tj"kxj^  -  d^»,k)J  dutdk  . 

Thus  we  find  a  set  of  simultaneous  equations  which 
must  be  solved  for  the  f  Clearly,  since  the  coefficients 
multiplying  f ^  depend  upon  (t^-t^)  and  (x^-Xj)  they 
have  the  same  symmetry  properties  that  were  stated  for 
multi-dimensional  autocorrelations.  If  follows  then  that 
once  the  integrals  have  been  evaluated,  the  recursive 
computation  procedure  given  in  Chapter  5.121  can  be  used 
to  solve  the  simultaneous  equations. 


Figure  6.2-3  shows  the  w-k  transform  of  three 
filters  competed  by  the  least-squares  method.  In  each  case 
the  weighting  function,  W(uu,k),  was  zero  everywhere  except 
in  the  PASS  and  REJECT  regions.  It  was  held  constant  ever 
each  of  these  regions.  In  the  three  examples  shown,  the 
ratio  of  the  weighting  function  in  the  REJECT  region  to 
that  in  the  PASS  region  (which  is  labelled  N/S  in  the 
figure)  was  set  at  1.,  10.,  and  100.,  respectively.  The 
remarkable  improvement  in  the  filter  performance  as  the 
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6.2 


Figure  o.2  -  4:  Coefficients  for  third  filter  (N/$  =  100.) 
illustrated  in  Fig,  6.2-3. 


ratio  increases  reflects  the  fact  that  we  are  much  more 
interested  in  the  reject  region  having  very  nearly  zero 
amplitude  at  all  points  while  we  are  less  concerned  with 
small  fluctuations  in  the  pass-region.  The  filtercoef- 
ficients  for  the  last  filter  are  given  in  Figure  6.2-4. 

As  a  final  illustration  of  the  effectiveness  of 
these  filters  when  applied  to  data,  a  set  of  traces  was 
constructed  which  contained  noise  whose  spectrum  is  pri¬ 
marily  in  the  REJECT  region  and  signals  that  are  in  the 
PASS  region.  These  traces  were  convolved  with  the  trans¬ 
form  filter  illustrated  in  Figure  6.2-1  and  with  the  least- 
squares  filter  illustrated  at  the  bottom  of  Figure  6.2-3. 
The  results  are  shown  in  Figure  6.2-5.  The  signals  had 
Btepouts  of  1.5,  1.0,  and  0.5  digitization  units  per  trace. 
Thus  the  signals  that  were  on  the  edges  of  the  pass  band 
of  the  transform  filters  were  attenuated  by  ~2.  On  the 
other  hand,  because  of  the  difficulty  of  constructing 
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Sample  of  noise  +  signal  traces. 


Noise  +  signal  traces  convolved  with  band-pass  filter 
(see  Fig.  0.2-2). 


Noise  +  signal  traces  convolved  with  band-pass,  band- 
reject  filter  (see  Fig.  6.2-4). 


Figure  6.2  -  5:  Example  of  application  of  band- pass  and 
band-pass,  band-reject  filters  to  simulated  noise  and 
signal  traces. 

noise  that  is  completely  restricted  to  the  REJECT  region, 
some  very  low  frequency  waves  were  actually  amplified  by 
the  least  squares  filter. 


167 


7 


utrtwiw  amWWlHI  IIIIWlii!;.'". . .  MlMMWi. 


APPENDIX 


PROGRAM  LISTINGS 

Nearly  all  of  the  subroutines  used  for  the  com¬ 
putations  made  for  this  thesis  were  performed  by  programs 
from  the  set  documented  by  Simpson  (1965).  A  complete 
description  of  the  writing  format,  the  abbreviations,  and 
other  programs  referred  to  will  be  found  in  the  work  cited 
above . 


All  of  the  subroutines  listed  here  were  written 
for  the  purpose  o**  investigating  polynomial  matrices  and 
their  spectral  decomposition.  Several  of  the  programs 
(BRAINY,  MATML3,  and  SIMEQC)  are,  however,  much  more 
general  in  their  application. 


nn  no  onn  mnnoonnnnnooonn  nnnonnonnnnn  on  on  onooonooooon  on  o  o  no  o  n  noon 


•RAINY 


Mncmc  LISTINGS 


BRAINY 


•  BRAINY  I SURROUT INC  I 

■  LRRCL 

CBAAINV 

SURROUT INE  BAASNV  I NRA ,NCA, l RA , * ,NCB , l OR  ,  »  ,C  I 
- ABSTRACT - 


TITLE  *  BRAINY 

COMPLETE  TRANSIINT  CONVOLUTION  Of  1*0  BATA  I  *  VAluEO  SERIES 

SUBROUTINE  BRAINY  COBRUTES  THE  COMPLETE  TNANSIENT 
CONVOLUTION 


Cl  II  •  SUB  OVER  All  VALUES  I  A I J- I » I  I •> I  I >  ) 

MHfRE  a,  B,  ANO  C  arc  BATRI  X-VALUfO  VECTORS. 

LANCUAGE  -  EUR  TRAN  II  SUBROUTINE 
EQUIRNENT  -  TON  OR  TONO  IBAIN  FRAME  ONLYI 
STORAGE  -  /AS  REGISTERS 
SREEO 

AUTHOR  •  R. A .  RIGGINS  ANO  J.E.  CLAERBOUT  |/A4 

- USAGE - 

TRANSEE*  VECTOR  CONTAINS  ROUTINES  -  SETRS. STf 
ANO  EORTRAN  SYSTEB  ROUTINES  -  NONE 

TORT RAN  USAGE 

CALL  8RAINVlNRA.NCA.lRA.  A.  MCfl.lBS.lt. Cl 
INRUTS 

NRA  NUMBER  ROMS  IN  A  MATRICES. 

NCR  NUBRER  COLUMNS  IN  A  MATRICES.  NUMBER  RONS  IN  B  MATRICES. 

I  BA  number  of  MATRICES  in  a. 

A I  1 1  1*1..  .NRA ,  I ...  NCR.  I...LRA  IS  A  CLCSElV  SBACEO  VECTOR  Of 

MATRICES. 

NCR  NUMBER  COLUMNS  IN  B  MATRICES. 

IBB  NUMBER  OE  MATRICES  IN  B. 

Bill  l•I...NCA,I...NC8.  I...LRB  IS  A  CLCSElV  SBACEO  VECTOR  OE 

MATRICES. 

•NOTE”  IE  NRA.  NCR.  NCR,  IRA,  UR  IRS  ARE  ISTHN  I.  BRAINY  RETURNS  Ml Th 
NO  COMRUTAT IONS. 


i 


CUTRUTS 

(III  l•I...NAA«l. .  .NCB.  I...LRA*'  BB-I  IS  TMC  CLOSELY  SBACEO 

VECTOR  OE  MATRICES  OF  THE  CDNVOIUTION  0*  A  Ml TH  R. 

EXAMPLES 

I.  INPUTS  •  NRA  •  t  NCR  ■  l  IPA  •  I  NCR  •  >  IBR  •  / 

All. .1.1. ..IT  *  l.t/., 4, 

Al I  ...  I . 1.. ./< I. . ./>  •  -  I • » • I. »- /  . .*4. 

OU  T  Ru  TS  *  Cll.../, I...?. I*. .4)  •  -  I • .-/• »- S ..-A. «-4. V. ,-|0. ,-l T. , 

-  T«  .*10.  i*  I  T  •  ,— / 4.  .-A.  .-B.  .  —  1 .  —  1 6. 


RR0GA4M  EOLLOMS  HEltl* 

OIMENSION  4ITI.BI2I.CIEI 
DIMENSION  s»acei/oi 

CALL  SETRS  I  NR  A.  M,  NCA.l  ,NCB.  M,  l  PA.l  A.l  PB.l  <1 
’  IXMlNOf  IN.M.l.lA.LBM  IOC, 100, S 
S  NT  I NUE 

CALL  SETRS  IM*N,MN,M«L,ML,L«N,LN,l,M|t|tR)l 
CALL  ST/  I  MN» 1 1 A*L B- 1 1 .C  I 
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•  BRUNT 

• 

MilCRRN  l  1ST INGS 

•  •• 

• 

BRUNT 

IRiCI  / 1 

IRRGt 

no  50 

L  1  -Ml 

00  *0 

*l«l  I 
K?«l? 

00  10 
JI«KI 
j)>K  I 

00  ?0  IR-I.N 
J  f»*f 
C«*0. 

00  10  IS*  1,1 
CK*CK»*<  J?I*R< Jll 
Jt*J ?*" 

10  J1*JI*I 

CIJI l*C**C< J 1 » 
to  Jt*JI*« 

«»•*! *1 
10  «?*R?*1 

H'lltWI 
*0  l?*l?**«l 

Ml -Ml  *«*N 
SO  M*M*IN 
100  Rf  TURN 
{NO 
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•  DIIMC 


MOSUN  LISTINGS 


•  •••••••••<>••••••••••••• 

•  OtlMtC  • 


•  CIHtC  (SUM  OUT  INC  I 


•  LABEL 
COI (MAC 

SUBROUTINE  OIINAC  (  U ,  1 1  ,NTR  A,  NCA.  AA .  *  I  ,BA  it  1 1 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


- ABSTRACT - 


TITLE  -  01  MAC 

CONnM  0 1  AGONAL  NATAI*  TINES  CONNIE!  SLUARE  NATRII. 


OilNAC  CONFUTES  THf  NATRII  NRQOUCT 


01  AC  I  l  I  •  A  ■  t 


NHERE  l  IS  A  CONNIE!  VECTOR  ON  LENGTH  NfRA. 

A  IS  A  CONNLEX  NIAA  BV  NC A  NATRII,  AND 
B  IS  A  CONNtti  n/RA  BV  NCA  NATRII. 


LANGUAGE  -  FORTRAN  II  SUBROUTINE 

EQUINNENT  -  TON,  TONO.  TON*  I Nfi |N  FRAME  ONLVi 

STORAGE  -  N)  REGISTERS 

SNEEO 

AUTHOR  -  R.A.  HIGGINS  I/A* 


-—USAGE - 


TRANSFER  VECTOR  CONTAINS  ROUTINES  -  NOT  ANV 
AND  FORTRAN  SVSTEn  routines  -  NOT  ANV 


FORTRAN  USAGE 

CALL  01  IN  AC  <!A./I.NlAA,NCA,AR,A!'BR'BM 


INNUTS 

INI  1 1  !•!,.. ..NINA  CONTAINS  THF  REAL  NART  UF  THE  VEC  >y0».  I, 

11(11  . . NINA  CONTAINS  THE  INAGINARV  NART  OF  THf  VECTOR 

I. 

NEAR  NUMBER  CE  ELEMENTS  !N  l  ANU  NUMBER  OF  ROMS  IN  A  ANO  B, 

NUST  RE  GRTHN>  I 

NCR  NUMBER  OF  COLUMNS  IN  A  ANO  B. 

MUST  RE  GRTMK*  | 

ARID  I  •  I  • .  . .  »NE  A  A»Nv  A  CONTAINS  THE  REAL  NART  OF  THE  NATRII  A 

STORED  CtOStlv  SPACED  BV  COLUMNS. 

Aim  l«l....,NERA*NCA  CONTAINS  THE  INAGINARV  NART  OF  THE 

NATRII  A  STORED  ClOSELV  SRACEC  BV  COLUMNS. 


OUT NUTS 

BRI 1 1  l-l,...,*fRA>NCA  CONTAINS  THE  REAL  NART  OF  THE  NATRII  t 

STONEO  CLOSEST  SNACFO  BV  COLUMNS. 

Rill  I  . . NIAAMCA  C'INTtINS  THE  INAGINARV  NART  OF  THE 

NATRII  B  STORED  CLU.ELV  SNACEI  BV  COLUMNS. 


El ANNIES 

I.  INNUTS  •  /R I  I  I  ■ fll|t«|.  N/R  A. <  NCA.J 
AR||...?I«I.,I.  All  |... *»•*., 0. 
OUTNUTS  -  RR( l...?l>0.,6.  R 1 1  I  ...  2  I •  *. , i. 


171 


oiamac 


PROGRAM  t  ISTINGS 


OUMC 


'MM  ?l 


IMAGE  ?( 


C 

c  2.  inputs  -  mi i. i.  /tu...2i*o..o.  hi* a.?  nca»i 
C  Ml  I. ..2I*I.M.  Alt  I. ,.2I«2.« *. 

c  outputs  -  bii  .■». 

t 

c 

C  PROGRAM  FOLLOWS  BElO* 

c 

c 

c  oup"v  dimension 
c 

DIMENSION  /HI? I ,1 || 21.AH12I , A  I  I 21 , A«  I ?>,BI <  ? I 

c 

C  BRING  IN  SC  AIM  VARIABLES. 

C 

N-N/RA 

•■NCA 

C 

C  FOR"  USEFUL  COMBINATIONS 
C 

NN*N»N 

m«i 

c 

C  00  MULTIPLICATION 

c 

00  20  1 1  ■  I  •  N 
00  10  I  AR«  1 1  . N"  t  N 

BRI I AB I  *  M I I / I • AR I  I  AB I  - { 11  IT ) • A  1 1 tABI 
HIIIATII'MI  IZMAII  IARI»/ II  12  MARI  TAB) 

10  CONTINUE 
20  CONTINUE 
C 

C  ThAT*S  ALL 
C 

RETURN 

END 
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L  AVEC 


MOMM  LISTINGS 


l*  VIC 


§ 


•  LAYEC  ISUAAOuTINEI 

•  Ltltl 
CL  AVEC 

SUAAOUTINE  LAVEC  ILI  ,EA, El. NACA.lPA, AOJA, EIFU, UA.ull 


- AHS  HAC I 


TITLE  •  LAVEC 

LATENT  VICTIMS  EO*  *  POLTNOWI  Al  MATA|I 


LAVEC  FINDS  THE  LE  LATENT  VICTIMS  Ul II  (M  VIII  OF 
A  FQLVNQNIAL  HATAU  AIEI 


T 

AO JUGATE  I  A  I  111)  I  I  •  UIII  •  VIII  . 


GIVEN  THE  COMPLEX  LATENT  ECAOS  Till.  1*1.... .11  ANO  THE 
*AT A | X  COEFFICIENTS  OF  THE  AOJUCATE  OF  AIEI.  THE 
VFCTOAS  AAE  SCALEO  SO  THAT  THE  F|*ST  NON-EE AO  ELEMENT  IS 
EQUAL  TO  I  I., 0.1. 


LANGUAGE 

EOUIAMENT 

STOAAGE 

SPEED 

AUTHOA 


FOATAAN  II  SUEMOUT INE 

TON.  701 C,  TON*  IMAIN  FAAME  ONLVI 

}?•  AEGISTEAS 

A. A.  WIGGINS  10/4* 


-USAGE - 


TAANSFEA  VECTOA  CONTAINS  AOUTINES  -  CHOOSE . 1*1 VEV.MATAA, 
ANO  FOATAAN  SYSTEM  AOUTINES  -  NOT  ANY 


FOATAAN  usage 

call  IAVEC  ILE.EA.EI.NACA.LPA.AOJA.EIFU.UA.UII 


INPUTS 

LE  LENGTH  of  VECTOA  OF  LATENT  EEAOS 

INI  1 1  . . LE  IS  AEAL  PAAT  OF  THE  VFCTOA  OP  LATENT  EEAOS. 

E.  OF  AIEI. 

Elll  I  . . LE  IS  IMAoINAAV  PAAT  OP  THE  VECTOA  OP  LATENT 

EEAOS.  E.  OF  AIEI. 

SAC  A  NUMAEA  OF  AOMS  OA  COLUMNS  IN  I  HE  PATAICES  OF  AOJ*. 

MUSV  At  GATHN*  I 

LPA  LENGTH  OF  POLYNOMIALS  IN  AOJA,  I.E.  THE  NUMOEA  OF 

COEFFICIENT  »ATA|CES  IN  AOJA. 

MUST  HE  GATHN*  I 

AOJAIII  l*I.....NACA*NACA*LPA  CONTAINS  T HE  MATA  I X  COEFFICIENTS 
OF  THE  AOJUCATE  OF  AIEI. 

EIFU  *0.  IF  THE  UIII  VCCTOAS  AAE  CESIAEO. 

NOT  *0.  IF  THE  VIII  VECTOAS  AAE  CESIAEO. 

OUTPUTS 

UAIII  1*1 .....NACAALE  contains  the  aeai  paats  of  the  latent 

VECTOAS  UIJI,  J*I.....LE.  IF  E!FU*0.  CM  VIJI  IP 
E  ITU  NUT*  0. 

Ul 1 1 1  1*1 ,....NACA*LE  CONTAINS  THE  |MArINAAV  PAATS  OF  Tm£ 

LATENT  VECTOAS  UIJI,  J*I,...,LE  IF  EIFll*0.  UA  VIJI 
IF  EIFU  NOT*  0. 
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IAVEC 


PBUGAA"  I  I  SUNGS 


LAVEt 


t }  I  I  PAGE  71 

c 

C  (UMtM 

c 

C  I.  INPUTS  - 

c 

C  outputs  - 
c 
c 

c  ?.  inputs  - 
c 
c 
c 
c 

C  OUTPUTS  - 

c 

C 

c  i.  inputs  - 
c 
c 

C  OUTPUTS  - 
C 
C 

c 

c  a.  inputs  - 

C  OUTPUTS  - 

c 
c 
c 

C  PROGRAM  E  01  (.OHS  At  10* 

c 

01  MENS  I  ON  7RI7I.7II7I.AOJAI  7  I ,UR I  7  I . U 1 1  7 1 
MNRCA 
l*LPA 

nn*n*n 

LN*l*N 

CALL  CHOOSi  I7IEU,  ILXlN.l  .IN.  llNXL.LN.il 
CALL  NAIRA  IAOJA.NN.L  .AOJAI 
AIM 

00  70  11*1,17 
?AC*?RI||I 

?|M.?|t||l 

J?*l 

oo  ao  i?*i ,n 

Jj*  J? 

RJ*R  I 

00  70  IJ*I.N 

CALL  IRLVIV  IL  .AOJAI  JJI, ERE. 7IM, URIAH. Ul  IK  Jl  I 
A)*A)»I 

70  JJ*JSMLXIN 

CALL  NAXAR  IN.URIKI I.UMX, |U*XI 
IE  I  ABSEIU*X  l-I.E-OAl  10.  SO.  AS 
JO  call  naxar  in.liiai i.umx. iumxi 
IE  IABSE IUNX l-| .E-OAI  AO, AO. AS 
AO  J7*  J7HINXI 
AS  CONI INUE 
AA*AI 

00  SO  I  A* I  ,N 
URL  *U« I AA I 
Ult*UllAAl 

ABSU-URl*URl»Ull*Ull 
IE  (ABSU-I.f-OAI  SC.SO.SS 
SO  AA*AA* I 
SS  CONTINUE 
AS*AI*N-I 
00  AO  l  A*  A  A  ,  A  A 
URT *UR I ISI 
U I  * *U I  I ISI 

UR  I  I S I ■ I URT  *UR l *U II*U|LI/ARSU 
AO  Ull ISI *IUIT*URL-URT*UIL I /AOSU 
70  A|*A|*N 

CALL  MATRA  I AOJA.L.NN. AOJAI 
RETURN 
ENO 


L  7  *7  7  Rl I... 71*1. ,7.  711  I ...  71 *t. ,0.  7IEu*0. 

NRCA*  |  lPA»7  A0J7I  l...7l*l..t. 

URI l...7l*l.,l. 

Ull  I  ...  7 1  *0. ,0 . 

17*7  7RI I... 71*2. .7.  7 1 1  I. . .71 •  I  . .- 1 .  7IEU*0. 

NRC A*7  L  PA *7  AOJAI I...SIH I..  0.1  I-.A66A7,  -.111111 

10.,  t.l.  I  • AOCOO  -.JJJJJI 

note  that  aoja  is  storeo  as 

AOJAI  t««.R!*t»,0.,0. «!«,-» A 46 6 A  7,.A,-.I||||,-.JJJJJ 
URI 1 • • • A I *| • « -0. A , I . , -  0. A 
U I ,  I . . . A  I  *0 . ,  I .8,0., -1.8 

17*7  7RI I . . . 7 1 •- I . . 7.  7111. ..71 *0. ,0.  7IEu*0. 

NRC  A*  7  L  PA*?  AOJAI l...il*l l.,0.  I  l-.S.O.I 

10. . 1.1.  I  .S.l.l 

URI  I...AI*!.,-.  IJJS.O. . I. 

Ul I  I  ... A  I *0.,0.,0*,0. 


SAME  as  Example  ).  EXCEPT  7IEUM. 
URI I.. .At*  l.,0.,l.,J. 

Ull  I ...  A  I  *  0. • 0. , 0 • , 0. 
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Mt#i  | 


MRO&RAN  1 1  STINGS 


•  NATNU 


•  MAT Ml l  I SURMOUT |Nf • 

•  f- AM 

•  MATMl  | 

COUNT  TOO 
l«L  MATMl  I 

ENTRY  MATMl  T  1  N,  M,  l .  A  A ,  RR ,  TENB  TA  ,  CC  .  G  l  E  AOO  • 


- ASST* AC T- 


Tlfll  -  MAT  Ml | 

GENERAL  NAT  A  I  *  MUl  T  I  Ml  It  AT  ION 


MATMl)  MULT  I  Ml  IIS  AN  N  ST  M  MAT* 1 1,  A,  IT  AN  M  MV  l 
MAI*|I,  H,  TO  OtITAIN  AN  N  IT  L  MMUDtiC  T  MAT*  I  * ,  C. 


M 

I  I 

VIA) 
t  I 


I 

I 

SI  M 

I 

I 


l 

I  I 
N  1C) 
I  I 


A  IS  ASSUMEO  TO  SI  STOREO  Sr  COLIMNS.  S  MAT  BE  STOMD 
EITHER  SV  COLUMNS  0*  *r  ROMS.  C  Mill  ME  STOREO  MV 

cm  urns. 


LANGUAGE 

EQUIMMENT 

STORAGE 

SMIEO 

AUTHOR 


E  AM  SUSROUI INI  IVOR  TRAN  II  COMMAIIHIEI 
TOR,  TORO,  TORN  (MAIN  ERAME  ONlVI 
TR  REGISTERS 

ASUUf  I M»»M,|0)»N*IOI«l  ♦  I T  0  MACHINE  CVCIES  UN  THE 
TORO. 

K.A.  MlGGlNS  R/TT/6A 


- USAGE - 

TRANSEEM  VECTOR  CONTAINS  ROUTINES  -  INOT  ANVI 
ANO  E OR  TRAN  SYSTEM  ROUTINES  -  INOT  ANV) 


EORTRAN  USAGE 

(All  MATMl |  IN,M,l,AA,B8 , IENRT* ,CC ,GIM AOO I 
IMMUTS 


NUMBER  CE  ROMS  IN 
■UST  HE  GRTHN*  1 

A 

ANO 

C. 

N>  MRE*  CE  COLUMNS 
MUST  BE  G«ThN»  1 

IN 

*  * 

ROMS  IN 

B. 

NUMBER  CE  COLUMNS 
•UST  RE  GMThni  1 

IN 

R 

ANO  C. 

l»l,...,N»M  CONTAINS 

1  Hf 

MATRI* 

All, Jl 

J ■ I • , . M  STOREO  CLUSEL  V  SMACED  HV  COLUMNS. 

SSI  I  I  1*1 . .  CONTAINS  The  MATRII  RII.JI  l«I...M, 

J - I  ...  I  STOREO  EITHER  8 V  RUMS  (R  COLUMNS. 

TENBTR  .0.  IE  RR  !S  STORED  RV  COIUMNS. 

NUT  *0.  IF  RR  IS  STOREO  MV  ROaS. 

GTE  AOO  GRTHN  0.  IE  IhE  MROOUCT  IS  TO  BE  ADDED  INTO  Ts*  OUT#UT 
AREA. 

ISTHN.O.  IE  IME  MRCOUCT  REPLACES  |H|  OUTPUT  AREA. 


OUT*UTS 

CCIII  I • I , . . . ,N*L  CONTAINS  THE  MROOUCT  MATRI*  CII.JI, 
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SUM'S.  1 


PP  )&■•  »“  t  1ST  (  Mi. s 


•UtNl  ) 


IP*  “  E  «M 


|0*f.»  /  S 


I  *  I • . *n.  stoato  ct:‘,*i*  sp»r*  «* 

Cl.'t  U“NS  . 


t  U»PUi 


1  NPSI?  S 

N  -  1  »»*  1  I’l 

A  A  I  1  *  >  . 

( 

Il*».  MNaiu  O.  GIUDU.O. 

outputs 

-  tuil'ft. 

*  N»U?  S 

-  **-  i  *•<>  t  -i 

A A (  |  .  .  .hi 

*  i 

. .  1 s. ,  1 . 

l  NRTK-Q. 

HIM  1  .  • 

*  i 

..S.. ). .1. 

Off  ADO* 1 . 

cr i i . . .6t 

•<  i 

. • 0. , 0. ,0. ,P. .0. 

Ob! PU  *  ^ 

•  CUI...6I  * 

1  t .  .  Ife. , R. 

.ii. 

,•>?.  ,16. 

INPUTS 

-  SA*«l  AS  i.  !«t»P! 

/ f \r t h  i.  r./raon-o. 

outputs 

-  Cf  1  I...M  * 

6.. 

)  . 

s*. ,??  . 

PAll'.PA*  IOllO»S  ‘ill  OH 


an  1 

HIP 

0 

*p/ 

hi  a 

c 

ads 

HIM 

0 

BCi 

1  .HAT  ML  ) 

sujml  j 

S*l> 

««*,N 

s«o 

»a  1 . 1 

s*n 

«*?.. 

Cl* 

s.s 

*00 

A  1 

St* 

** 

Cl  * 

S.s 

»D0 

K  1 

si* 

R« 

C!  * 

I.S 

*nn 

A  1 

st» 

IPS 

St* 

cc 

St* 

CCI 

IOC 

NOP 

Cl  *  • 

S.s 

CAS 

/IPO 

•  P* 

•  •1 

NOP 

NOP 

log 

St/ 

SlO 

t  pp 

Cl*» 

1 .6 

SlO 

u 

sto 

Nl 

*c» 

UP*. 

/ » * 

**S 

l 

Ps« 

.  1 

S«* 

ent .  1 

Cl  A* 

iog 

nu 

/FT* 

fc.S 

*r.« 

SID 

hupo 

»CA 

SID 

OOP  N 

100* 

/.  * 

UP*. 

l.s 

*1  s 

1 1 

sto 

pot 

iog* 

: .  s 

■  P*  • 

».* 

*»s 

i 

p*« 

.* 

Nl 

1 1  * 

•  *  / .  i .  • 

All 

0.  1 

S«U 

N"l  .  1 

l  « 

PCI  .  1 

i  so 

Pit./ 

m* 

*  PH 

ICOP 

<  1  * 

N**|  .  1  .  1 
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F«0&*  »«  C  I  St  INC, S 


i 


•  F*f«Lt 

jp«ce  >i 


X* t«l i 


IFFGf  1) 


M* 

eot  .i 

MO 

80t  ,7 

■0*0 

t  l« 

1  ,7 . 

•  <)*  ! 

S'O 

not  ,7 

Ttt« 

l  *■ 

N«l 

f  II 

•  ••*•<•*  l) 

l  <0 

8  It  ,7 

IF# 

S  f  7 

•*.* 

•••A0H1CC) 

It 

L  JQ 

1 

•  ••AC*  < HA  | 

HA 

M«F 

••.7 

•  •■AC* 1  A*  » 

CC 

MO 

•  •.AO«<CC  » 

CCI 

stu 

•  •»  AD*H  CC  1 

ao«« 

in 

I 

•••I  Q*  * 

N 

t  ii 

Hi  1  .  »• 

•  • 

tu 

t  OOF,  * ,  | 

MO 

*«l,l 

L  >0 

*•7.7 

MD 

ll«,l 

tft* 

HOI 

F7f 

••■N*M, *  •• 

/MO 

F7C 

Sf  7 

S  T  7 

•  • 

KID 

F/f 

..1 

Kl 

Ptf 

1  ,  • 

f  NO 
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1 


POhAIN 


PAOGA AH  t  (STINGS 


PQHAIN 


•  POHAIN  ISUMAOUlINfl 

•  L*»fl 
CPO**IN 

SUOAIXMlNt  POHAIN  INACA.LPA.A.i  40J.4UJUG.0l  TPOl  .SMCIl 
C 
C 
c 

c  mi*  -  domain 


- AASIAAC- 


PniYNOHIAl  HAIA|A  AUJUGA I f  AND  UftfMH|N«Nt 

sosaoutinc  po«ain  f  inos  iHf  »ojoc*if  and  otitAHjNANi  uf 

A  l AN ADA- HAT A  1 1  IA  L AHPOA- HA  I A  I  I  IS  HAIAIA  IN  hh|Cm  f ACH 
IfAH  is  A  POLYNOMIAL,  OA  fQUIVALtMlv,  II  ISA  POLYNOMIAL 
HAVING  HAInII  Cliff  fit  If  NISI  .  IHf  Nf  TmOO  USIU  IS  AN 
f  A  I  f  NS  I  ON  n*  IHf  I  NVf  A  S  ION  ICCHNUHf  UfSCXIAfO  e» 
f  A  OUt  Jf  H  AND  SON  |NS»  |  I  PXOBLfH-AOi  «  IN  h  |  r.Hf  A  AlGtBHA, 
SfCONO  f  D  1 1  I  ON .  GQSIfCHISDAT  I  <f  A  V  I  .  THIS  MTfNSIUN  HAS 
HAD!  At  J.  CL  A  f A  SOU  I .  STAIIS'lCS  INSIIIult.  UNlVfASItY 
Of  OAASAIA,  SHfOCN.  AND  A.  JANS  SO  \  •  Af  St  A  A  f  H  INSIHUff 
Of  NAllUNAL  OfffNSf,  ifOCAHOLN,  SkfDfN. 

-  fOAIAAN  II  SUrtNOUflNt 

-  FDA  OA  FOAO  I  HA  IN  fAAHf  ONI  VI 

-  212  ACGISItAS 


C 
C 
C 
C 
C 
C 
C 

c 

C 

c 
c 
c 
c 

C  LANGUAGC 
c  (quipmint 
C  STOAAGf 
C  SfffO 
C  Au  I MCA 

c 
c 
c 

C  IA.  NSf  f  A  VtCTOA  CONTAINS  AOUllNtS  -  BA  A  IN  V  ,  HO  V  C  .  SI  I A  S 


-  A. A. 


HIGGINS  ANU  BO  JAMSSOM  i/** 

- USAGf - 


ANO  fOAIAAN  SVSTfN  AOUllNtS  -  NONf 


C 
C 

C  fOAIAAN  USAGf 


C 
C 

C  INAuIS 
C 

C  NACA 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C  OUlAUfS 

c 


CALL  AnNAlNINACA.i»A,A,LACJ. AOJUG.Dt lAOl.SAACf I 


l  PA 


All' 


NUNHH  Of  A  OHS  OA  COLOHNS  Ml  IHf  PATAICfS  Uf  A  . 

HuSI  Pt  G4lHN«  I 

LfNGtH  Of  I  Mf  POLYNOMIALS  IN  A  I  I  .f  .  IHf  NUHtfA  Of  HA  I A  |  A 
POLYNOMIAL  COtf  f  1C  If  NIS  IN  Al. 

HuSI  At  GAlHN«  I 

!•!... NACA. I...NACA. I...LPA  |S  IHf  VtClOA  Of  HA  I A | A 
CntfflCIfNlS  f OA  IHf  POLYNOHIAl. 

HUSI  HI  SIOAfO  CLOSfLY  PACAfO. 


SPACt  I  I  I 


I >1 . . •  N  AC A»NAC A • | I NACA I. 1 1 PA- I  I* l  I 
NffOCO  PY  POHAIN. 


is  roNPuiAi ion  sPAtf 


LADJ 


AD JUG  I  1  I 


c 
c 
c 
c 
c 

t 

c 
c 
c 

C  f  A AhPl  t  S 

c 


LfNGlH  of  POLYNOMIALS  in  AiJJUGATl  Of  A. 

•  I  NACA  HUPA-  I  I  •  I 

l» I ...NACA.  I., .HACA.I...LAUJ  IS  >f  VtCTOA  Of  HA  I  A  I  A 
COffflCIfNlS  f  0»  IHf  AO  JUG  A  I  f  Ok  A. 


OffPOtlll  l»l ...LAOJUPA-I  IS  IHf  VtCTOA  0»  COffflCIfNlS  Of  IHf 
P0LYN(,“|Al  OflfAHlNANl  Of  A, 


C  I. 
C 
C 
t 

t  2. 

C 

c 
c 
c 
c 
c 


I NPU I S 
OUIPUI S 


INPUIS 

OUTPUTS 


NACA  •  | 

-  LADJ  •  I 
Dt  I  Pill  *  I  . 

-  NACA  •  2 
Al|...f.l 

-  LADJ  •  2 


L  PA  >  2  All. ..I, l, 

ACJUGI 1 • • • 1 . I ...I. I. ..I  I  *  I. 
.  .2 1  •  I..I. 


I..I. 


I  Pa  •  2 

I...I)  *  l«.li.r..la|L).L..L.il. 
Of  I  POL  II...  II  •  I..I..IJ. 


AP  JUG  I  I «.«/. I. • • I . . *  2  I 


.-I 


.  I . 


i.  INPUTS 


-  NAC  A  •  2 

Al  I  ..  .2,  I 


L  PA  •  1 

*.?.l«.«ll 
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ponain 


PAUGKAH  i  I  S T | NGS 


ponajn 


IPAG*  i  I 


I  PAGE 


c  outputs  -  l AO j  •  i  oetpoliii  • 

C  AOjuGl  l,«,?«l,,a?»l,,,II  * 

c 

C  PAOGAAN  fOLLU«S  BELO* 

c 

01  ME  NS  ION  M  8 1 , AO JUG  I  8  > ,  Of TPOl I ?) , SPACE  I  81 

CALL  Sf  TaS  I NACA.N.lPA.LAI 

CALL  SHAS  IN»N,NN,NM,NI,LA,lADI 

CALL  SETaS  1  NN  •LAvNNLAfNNLA,  NNl  A  1  ) 

If  IN-11  5,5,6 

5  CONTINUE 
l  AO  J»  I 

AO  JUG  I  I  I  •  |  . 

CALL  NOVf  I  L  A,  A i Of  T  POL  I 
GO  TO  80 

6  CALL  NOVf  l  NNl A , A, SPACE  I 
?000  1*1. 

DO  50  L • I . 5 

C  CALCULATE  COE  ff  I C l ENT  S  P1AI  Of  C"AA  AC  TE  A  1  ST  1C  POLYNOMIAL 
KIM 

CO  i 0  KM,NNLA|,NN 
PlA-O. 

NNl »A»NNM 
00  10  I  «K ,  NN 1 ,  N I 
10  PlA«PlA.SPACCI  1  I 

Cf TPOl I  A  I  I .PlA/f LQATf 1 1  I 
}  0  A  1 »A  I  M 

If  ll-Nl  ?5 , 50, 50 
}  5  CONTINUE 

C  SUBTMACT  PIKIMOENTITY  MAMMA 

CALL  move  (NNIAI  , SPACE. AOJUGI 

00  >0  I  »  I ,  NN ,  N  i 

AIM 

00  JO  AM,NNlAl,NN 
AOJUGI A  I • AOJUGI A l-OETPOl I«1 I 
JO  A I  *  A I  *  | 

C  multiply  by  input  matai* 

CALL  BAA INY  I N, N,l A, A ,N, L A£, AO JUG , SP ACE  I 

JOIO  ?•!. 

LAO’L  AO*L  A- I 

50  NNl  A|  >NNLA1*NNIA-NN 

C  CHANGE  SIGN  Of  OETEAMINANT  ANO  AO  JUGATE  If  N  *  VC  N 
50  CONTINUE 

If (IMOOf IN.JII  60,60,70 

60  CONTINUE 

00  61  IM ,LAO 
6|  OETPOLI  I  I  •-OCT  POL  I  II 
00  6?  I • I, NNL A  I 

6 1  AOJUGI I (••AOJUGI  I  I 
TO  CONTINUE 

IAOJM  AO-l  A*  | 

80  AETUMN 
ENO 


i  I 
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ti|)(  IttiHijiiirtiiiiiiWIItliMiiHiiiiwiiwc-MiM."  . .  . . .  1llll|lltoiiiWlilll|!il|J»:«®||!|||!||||||||j|t5ili(|l 


PROGR  AM  L  |  SUNGS 


SI"EQC 


St"f  OC 


SINEQC  (SUBROUTINE! 

LABE  l 

>  I  "E  QC 

SUBROUTINE  SIMEOC  < NR AO  I M, NR AB ,NCB , AR , A  1 , RR , R | , 01  T , I S . I R A  1 


TITLE  -  SPJEQC 

SOLUTION  OE  COMPLEX  SIMULTANEOUS  E3JATHNS 


SIMEOC  SOLVES  T ME  CUMPLEX  SIMuifXNEOUS  EQUATIONS 


AX  •  H 


c 

FOa  I,  KMfHf  A 

HAS 

Aft 

ROMS 

AND 

NRAB 

COLUMNS. 

c 

ft 

HAS 

Nil  A  B 

ROMS 

ANO 

NCB 

COLUMNS, 

c 

X 

HAS 

N«  Aft 

ROMS 

ANO 

NCR 

COLUMNS. 

The  SUL  UT I  ON  MATRIX,  X,  IS  STORED  IN  A. 

THE  SOLUTION  OE  THE  MATRIX  EOUA T I IN  IS  ACCOMPLISHED  BV 
UPPER  TR I ANGUL AK | / AT  TON  OF  THE  A  MATRIX  USING  A 
MODIEIEC  MAXIMUM  PIVOT  IONLY  THE  INDIVIDUAL  COLUMNS  ARE 
SEARCHED  FOB  A  MAXIMUM)  FOR  EACH  REDUCTION.  THE 
DETERMINANT  i  S  COMPtlTEO  AT  THT  SAME  TIME.  IF  T  hE  MATRIX 
A  IS  SINGULAR  A  TER U  VALUE  0»  Tml  DETERMINANT  IS 
RE  LURNIC 

BOTH  THE  A  AND  H  matrices  are  DESTAOVED  By  SIMEOC. 


LANGUAGE  -  FORTRAN  II  SUBROUTINE 

EQUIPMENT  -  TOR,  TORO.  TDRM  (MAIN  FRAME  ONL V I 

STORAGE  -  ATS  REGISTERS 

SPEED  -  RfLATiVEl  *  SLCJM 

AUTHOR  -  «.A.  mIGGINS  TTSM 


- USAGE - 


TRANSFER  VALIUM  CONTAINS  ROUTINES  -  NOT  ANV 
ANO  FORTRAN  SYSTEM  ROUTINES  -  NOT  ANY 


FORTRAN  USAGE 

CALL  S IMF  QC  INRAO|m.nRAB.NCB,A»,A| ,BR,B| ,01 T , IS.ERB? 


INPUTS 

NRAIMM  THE  DIMENSION*  0  SI/F  OF  THE  COLUMNS  OF  *HE  A  AND  B 

MATRICES.  That  is,  the  CALLING  PROGRAM  moulo  contain 
statement*  lire 

OIMfNSION  AR I SRAOIM, IGNORO) ,  A| INRADI m, ignoroi 
DIMENSION  RAINRRDI", IGNORCI ,  B| INRADI". IGNOROI 
"HERE  IGNORO  may  re  ANY  CONSISTENT  VALUE.  IF  TMt 
MATRICES  ARE  STUREO  CLOSELY  SPACED  THEN  NRA0I"  •  NR LH 
ANO  AR,  Al,  BR,  AND  B|  MAY  PE  DIMENSIONED  AS 
VECTORS. 

MUST  RE  GRTHN«  NP.AA. 

NRAB  number  CF  ROMS  IN  A  ANO  B,  NUMBER  OF  COLUMNS  IN  A. 

MUST  ftF  OR  T  HN*  I  . 

NCt*  NUMBER  OF  COLUMNS  IN  t) . 

MUST  BE  GATuN*  I  . 

ARII.Jl  1*1 . NRAB,  J • I , .. . ,NRAH  CONTAINS  THE  BEAL  PART  OF  A 

STORED  BY  COLUMNS. 

NOTE-  AR  IS  ALSO  AN  OUTPUT. 


A I  I  I  Jl  1*1, ....NRAB.  J*L,...,NRA8  CONTAINS  I  HE  IMAGINARY  PART 

OF  A  STORED  BY  COI  UNNS 
NOTE-  Al  IS  ALSO  AN  OUTPUT, 
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. . . ••••••••  Mocun  iisthCs  . . . 

•  SIMRQC  •  •  SIMRQC  • 

IMAGE  2  I  ( PAGE  2 1 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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AAll.JI  l*l,...,NAAB,  JM.....NCB  CONTAINS  ThR  REAL  RART  OR  » 

STORCO  AT  COLUMNS. 

IS  OCTROYED  BY  St"E<JC. 

BlII.Jl  . . NAAB,  . ,NCB  CONTAINS  THE  IMAGINARY  RANT 

OR  B  STORCO  BY  COLUMN* 

IS  DESTROYED  BY  SMCOC. 

SRACE III  Ml,. ...NAAB  IS  TEMPORARY  COMMUTATION  SPACE. 


OUTPUTS 

ARM,  Jl 


. . NRAB,  J*  I . NCR  CONTAINS  THE  REAL  RART  OF  R 

STOREO  BY  COLUMNS. 


AIII.JI  1*1  ,...,NRAB,  . . NCR  CONTAINS  THE  IMAGINARY  RART 

OR  I  STORCO  BY  COLUMNS. 


OEM  I  I 


1*1 . 2  CONTAINS  THE  REAL  ANO  IMAGINARY  PARTS,  AISPCC- 

TIVELV,  OR  THE  OETERMINANT  OR  A. 


ERR 


•0.  IF  ALL  OR. 

•  2.  IF  A  IS  SINGULAR  10CTI1...2I  •  0.1. 


EIAMRICS 

1.  INPUTS  -  NAAOIN.?  NAAB* l  NCB*2 

AMI  ll-2.  AllH*l.  BR 1 1 .  . .21 •/.  .  ]«  Bill. ..21"'  .-2. 
OUTPUTS  -  ARI I.. .21-1.. -B  A 1 1 1 .. . 2 1  >0. , - 1 . 4  OR  T I t. • .21*2. . > » 
ERR-O. 

2.  INPUTS  -  NRA0IM.2  NRAB-2  NCB*2 

ARI  I  •••A|m|..*I.,*I.,*I.  BR ll..«4l*I.,0. jO.,1. 

A 1 1  I . • .41 • | • ,  I*.  I. .-2.  B 1 1  | . . • 41 »0. ,0# , 0* .0. 

THAI  IS,  A  •  I  I  t.,  1.1  I - 1 . .  1.1  I 
I  1-1..  1.1  1-1..  2.1  I 

OUTPUTS  -  ARI |... 41— 0.1,  2..  I..  0.  CETII...2I  •  I-I..I.I 

All  |... 41— 2. 1,-1.,  0.,  I.  ERR  *0. 

THAT  IS,  *  •  I  I-. 1.-2. 11  II. .0.1  I 
I  l  2. .-1.01  10., I. I  I 


1.  SAME  NUMBERS  AS  MUMPER  2.  BUT  M|Th  A  OIRRRRRNT  0 1  ME  NS  ION. 

INPUTS  -  NRAOlMM  NRAB*  2  NC8-2 

ARI I...AI*!.,— I.,0.,— I.,— l.,0. 

All I...Al*l.«  1 . 1 0 . ,  I.,  —  2..0. 

BR  I  I*. .AIM.,  0  .  ,  0  .  ,  0.  «  l.,0. 

All  I  ...  A  1*0 . •  0..0.,  0.,  C.,0. 

OUTPUTS  ARI I...AI— 0.1,  2. ,0. , I . , 0. . 0.  CETII...2I  ■  l-t.,1.1 

All I...AI— 2.1,l.,n.,G.,l.,0.  IRR*0. 

4.  SECONO  ROM  IS  COMPLEl  CONJUGATE  OR  FIRST  AO*. 

INPUTS  -  NAAOIN*?  NR A9*2  NCB* I 

ARI I.. .41*  1.,  I.,  —  I.«*l.  BR 1 l . . .2 1 •  I.,  2. 

Allt. ..41*  |.,-l.,  l.,-l.  Bill. ..21*  1..-2- 

THAT  IS,  A  ■  I  II..  t.l  1-1.,  I. I  I  B  ■  I  II.,  2.1  I 
I  II. .-1.1  i-l.,-1.1  I  I  II., -2. 1  I 

OUTPUTS  -  ARII. ..21*1.1,0.  OCT  •  10., -4. I  ERR*0. 

All  1 . . .2  I  ■  0...S 


1.  SINGULAR  CASE. 

INPUTS  -  NRA0IM.2  NRAB»2  NC8*I 

AR 11. ..41*  I., -I.  -Mil.. .21*  I  .  ,  0  . 

AIM. ..41*  l.,l.,-|.,-L.  Bill. .*21*  0.,0. 

OUTPUTS  -  OR T I  1.. .21*0. .0.  E«R*2.  AR.At  CCNIAIN  MEANINGLESS 

NUMBERS. 
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r>  o  o  o  o  o  o  o  r»  r* 


SINtOC 


PAO&AAP  l  I  STINGS 


si  mot 


IMG*  tl 


IMG* 


C 

01  Nf  NS  I  ON  AA  I  7  I,  All  7>.BAI7l  ,811  71  ,0*  Tl  II,  I  Slit 
C 

C  MING  IN  SCAlAA  v*«  |  ABLE  S 
C 

NO*N<UOIH 

N*NAAB 

A*NCB 

EAA*0. 

C 

C  I*  NAA8  •  I  COMPUTE  OIAECUT  ANO  LEAVE 
C 

If  IN-(!  7. 7, A 
7  CONTINUl 

Otllll>M 
Clf  I2IUI 

Of  T  M*  I  Of  T  II  I  *0E  T  I  1 1  *0t  T  <  7 1  *DE  T  <  7 1  I 
J*( 

00  )  I  •  I  ,  A 

AAfJI*(AA(JI*OfT(ll*BIIJI*OETI7ll/OffM 
A  I I JI*IB| | J) #0* Tl I  I -BUI J I *0* 1(21  I /Of  T» 

T  J*  J*NO 

GO  TO  SOO 
A  CONTINtlf 

COMPUTE  THE  NffUEO  COMBINATIONS  ANO  SET  UP  THf  INITIAL  VALUES. 

NNO*N*NO 
A  NO*  K • NO 
01 T I  11*1. 

oeti?i*o. 

I  TA  •  I 
NOItftl'O 
00  10  I  ■  I  »N 
10  ISIII'I 

I  T  A  *  I  NOE  A  Of  T  A  ACT  TEAM 
NO  I T  A | *N0*<  I T“  -  |  I 


F I  NO  LAAGEST  VALUE  IN  THIS  CULO«N  KITH  I  NOE «  GmThn*  I TA 

100  CONTINUE 

IA* I TA 

AHA  *0 . 

110  CONTINUE 

I  4  I  A  «  |  S  I  I  A  I 

IS  I  A  A* IS  I A*N0l  f  A l 

AM*  T  ■  AM  I  I  S  I  a  A  I  •  AM  t  |  s  I  A  A  I  *  A  |  I  |S(AAI*AI(  I  4  I  ■  A  | 

if  (AMi-APiTi  no.  110. 110 
no  continue 

AMI* AH«T 
ICC  » IS  IA 
ISI  I  A  I  * !  S I  I T  A  I 
ISIITAI* ICC 
no  CONTINUE 
I  A*  I A  *  I 

if  iia-ni  no. no, no 
no  continue 

ICC  -  I S  «  ITAI 
2 000  0*1. 
c 

C  0 1  V  !  Of  ALL  TEAMS  TO  THf  A  I GHT  Of  ICC  BT  UlAGONAl  TEAM 
C 

ICC  A* | CC  *NU 1  I A  I 
AA  |  TA*AA  I  ICC*) 

A |  | TA.AI I  ICC Al 

700  CONTINUE 

ICC  A* | CCA*ND 

If  (ICCA-NNOI  710,710.770 
710  CONTINUE 

AAICCA.AAI ICCAI 
A |  | C  C  A  *  A | I ICCAI 

aAIICCAI •( AA ICC A* AN  IT A* A | I  CCA* A II TA I /A«* 

AinccAi*m  :cca*aai  ta-aaicca*aii  t  a  i  /  a**  « 


)l 
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SIAEQC 


PAOGAAM  L  1  S  T  |  NO  S 


SINfQC 


IP  AO* 


91  IMAGE  A) 

GO  TO  >00 
210  CONTINUE 
ICCB-ICC 
2  ‘.0  CONT  I  NUE 

8AICCS*BM( ICCSI 
Ml t CCS* fl  M ICCSI 

8* I  I  CCS  I *ldAICC8*AAIfA»BI 1C CB* A 1  I  TA  I /AN* 

81  I ICC8I » ( A | | CCS* AH  I TA-88 ICCB« Al ITA I /AMI 
ICC6*  ICCA.NO 

it  iicca-kni>>  >io.>io,>ao 

>80  CONTINUE 
>010  C*l. 

INCAENENT  (HE  UE  TEAM  INANT 

Of  T  |  *OET  I  I  ( 

CE  T>«Of  T  1  >  I 

Of  T(  H«0EM»A8|T8-niT>»A|  ITA 
Cf Tl >1 «0U I »A| |TA»OET>*AA I TA 

SUAf  MAC T  THIS  A0«  TAOM  All  SUCCEEOING  MOMS 

I A» I TA 

>80  CONTINUE 

>0>0  Q>l. 

IA* I K • I 

IT  IIA-N1  >60. >60.  1(0 
>60  COM»l HUE 

I  CCA* IS! I A ( • NO i I A  I 
I CC  S  *  I CC  *N0 1 T A  l 
AA ICC* AA I  ICC  A ( 

A!  |  CC*  A|  (  ICCAI 
>  TO  CONTINUE 

ICCA.ICCA.no 

IT  IICCA-NNOI  >80, >80. >90 
>80  CONTINUE 

ICCS*1CCS*N0 

AA  I  ICCAI  *  AA  (  ICC  A  I  •  I  AM  (  !  CCS  ( •  AM  ICC- A 1 1  ICC  S I  •  A 1 1  CC  I 
All  ICCAI *A((  ICCAI- 1 AA|  IClC<*AllCC*AlllCtSI*AAICC> 

co  ro  > to 

>90  CONTINUE 

ICC8*  |  S I  I A  I 
(CCS.ICC 
100  continue 

AA |  ICCAI* AN  I  ICC  A  I -  I  BA | ICCSI*AAICC-B|( I CC  S 1 •* I IvCI 

A 1 1  I  CC  A 1  ■  B 1 1  ICC  A  I  - 1  BA  I  I  CCS  I  ••  1 1  CC  *A  |i  |  CC  SI  •  AM  ICC  I 

(CCa*ICCA*NO 

<CC$*ICC$*«0 

IT  I ICCB-A NO  I  100, 100, >80 
110  CONTINUF 

LOOP  TO  NEAT  AON 

|TA*|TA»| 

NO|TA|*NOITA|«NU 
IT  IITN-NI  100,  100,  1>0 
1>0  CONTINUE 

IT  OETEAminANT  IS  EEAO,  LEAVE 

IT  I AAST (OEI I | ll*AASTIOf Tl >1 l-l.f-6)  999,899,110 

110  CONTINUE 

non  ClEAA  The  UPPEM  (AlANGlE 
ITA.n 

NO  I  T  A  I  *NNO-NiJ 
190  CONTINUE 

I  A* | T  A 

16S  CONTINUE 

>010  C*l. 

IA- I  A- | 

IT  I  I A  I  170,1  TO,  180 
180  CONTINUE 

ICCB-ISt (A| 
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SI**ICC 


PM  IGRAM  list  INGb 


SIRTCC 


IR6« 


S)  I  PACT  M 

ICCS- I  SI  I  Tm I 
ICC*" ICCR*NOI IH| 

*R I cc ■ an i icc  41 

4  1 ICC-41 ( ICC  *> 

160  CUNTINmt 

RRI  I  r.  C  R  I  *RM1  ICCR1-I8RI  ICCS1*RRICC-RII  ICCSI  *41  ICC  1 

n  1 1  icc  n  •  •  n  1 1  iccri-i nil  iccs i*ar icc  *rri  t  ccsi  *a  i  icr  t 
ICCB* ICC*»NO 
ICCS*  ICC$*NI> 

11  1ICCS-RN0I  160.  160. 

110  CONTINUt 
ITR«|TR-| 

N0ITR|*N0|TR!-ND 
11  HTR-II  ISO,  180.1*0 
180  CONTINUt 

UNSCRA«Blf  8  INtU  * 

IR»  I 

NO  I C I *0 
140  CONTINUt 

I R'NDICI *1 
400  CONTINUt 

I S  I R  A*  I S  I  I «  1  <  NO  I C 1 
*RI  IAI*8»» ISIRAI 
*1  I  1*1 «RI ! IS  I R 4 1 
|  A  «  |  4  •  | 

!*• I M* I 

II  IIR  N1  400.400.410 
410  CONTINUt 
IR*  I 

no  i  c  i  *  '40 1  c  i  *no 

II  INOICI-RNOI  140. 410, 420 
420  CONTINUT 

THAT  '  S  All 


soo 

CONTINUT 

Rf TURN 

994 

CONTINUt 

IRR ■ 2 . 

GO  TO  400 

T  NO 

ld^4 


SNLMAC 


MOO RAN  LISTINGS 


SNLMAC 


•  SNL MAC  I  SUBROUTINE) 

•  LABEL 
C SMI  MAC 

SUBROUT  I  NE  SNL  MAC  I NZ , ZR , Z I , NRU , UR , U I , L* . A • . SPACE .ERR  I 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

p» 

m 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


—  AASTAAC  T - 


nut  -  SNLMAC 

synthesize  L*»80A  MAT* 1 1  FAOM  COMMIE*  VEC  TOMS  AND  ZEROS- 


SNL  MAC  CONSTRUCTS  *  LAMAOA  I  POL  Y KOMI AL I  MATAI1  FROH 
THE  COMPLEX  LATENT  VECTORS  Ul I  I  ANO  LATENT  ZEROS  Zlll. 
SINCE  THE  ROtYNOMIAL  MATRIX  MILL  rtf  REAL,  THE  NON-REAL 
l  AUNT  VECTORS  ANO  ZEROS  MUST  APPEAR  IN  CONJUGATE  PAIRS. 

IF  ME  FORM  THE  MOOAL  MATRIX 


I 


U 


N 

<U(I  I  I  I 

i urn  i  i 

.  I  M 
.  i 
I  U I  M  |  |  | 


ANO  THE  SPECTRAL  MATRIX 

Z  •  01  AG  I  Zll >,ZI?ll....Z|M|  | 

THEN  THE  N  8Y  N  MATRIX  COEFFICIENTS,  Alll,  . . . 

ARE  GIVEN  8V  The  SIMULTANEOUS  EQUATIONS 


L 

AI|I«U«I  ♦  ...♦  AILI*U*Z  •  U. 

MHCRE  L*N  •  M  . 

SNL "AC  SETS  UP  ANO  SOLVES  THESE  EQUATIONS. 


LANGUAGE  -  FORTRAN  II  SUBROUTINE 
EQUIPMENT  -  ION.  1090.  1094  INAIN  FRAME  ONI Y I 
STORAGE  -  198  REGISTERS 

SPEED 

AUTK"1*  -  R.A.  HIGGINS  8/A4 


—USAGE - 


TRANSFER  VECTOR  CONTAINS  ROUTINES  -  OIINAC,  NAIRA,  NOVREV,  SINEQC 
ANO  FORTRAN  SYSTEM  ROUTINES  -  NOT  ANY 


FORTRAN  USAGE 

CALL  SNIHAC  INZ.ZR.ZI.NRU.UR.UI.LA.AA.SPACE.ERR) 


INPUTS 

NZ  NUMBER  OF  LATENT  ZEROS  ANO  VECTORS. 

•  M  IN  THE  ABSTRACT. 

MUST  BE  GRTHN*  I 

• 

ZRIII  1-I.....NZ  CONTAINS  THE  REAL  PARIS  Of  THE  LATENT  ZEROS 

Zlll. 

■ 

Zllll  1*1.. ...NZ  CONTAINS  THE  IMAGINARY  PARIS  OF  THE  l ATEN’ 

ZEROS  Zlll. 

N*U  NUMBER  OF  ELEMENTS  IN  EACH  VECTOR  Will. 

MUST  BE  GRTHN*  I,  MUST  BE  A  MULTIPLE  Uf  NZ. 

• 

UR  1 1  I  l-l . .  CONTAINS  THE  REAL  PARTS  OF  THE  LATENT 

VECTORS  Ul  I  I  STilRFO  CLOSELY  SPACEO. 


| 

i 

3 

1 


3 
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. . . . . 


PROGRAM  l  |UhGi 


SHI MAC 


I  PAGE  ?1 


UIMI  I  'I  ,..t,NRU*Nf  CUNT*  I  NS  THE  IMAGINARY  PAR  T  S  OF  THE 

LATENT  VECTORS  Ufll  STORED  CLCSElY  SPACED 

•THE  LATENT  VECTOR  Ul  It  CORRESPONDS  fU  THE  LATENT  AOOT  fill. 
both  VECTORS  MUST  CUNTAtN  The  Cl'PPLM  COFJUOATf  ■>  MPlIClTir. 


SPACE  I  I  I 


I • I  .... • ? *N/ •? • INf »NRU1 *NRU  IS  TEMPORARY  COMPUTATION 
SPACE. 


I  PAGE  fl 


C  OUTPUTS 
C 


IS  the  NUMBER  OP  COEFFICIENT  matrices  found. 

•INRU/Nf 1*1 

*  1*1  IN  THE  AbSTAACT. 

1*1 . .  A  CONTAINS  TMf  COEFFICIENT  MATRICES 

Alii,  1*0. ..I,  STORED  CLOSELY  S»ACEO  BY  COLUMNS  HHERE 
AIOI  •  IDENTITY  MATRIX. 

*n.  IF  ALL  OK. 

•7.  IF  The  SIMULTANEOUS  EQUATIONS  ARE  SINGULAR  -  THIS 
HAPPENS  IF  III  I  I  *Ut  J  I  WHEN  flll*f|j|. 


C  EXAMPLES 
C 


-N I’J  fRI  I. ..?»•?. ,J.  f  II  l...7l*C.,0. 
NRU'I  URI  I. ..?.*!. ,1.  Ull I . . ,  71 *0. ,0. 


Outputs  -  la*  )  aai  i... »i *i .,-s..b.  erp«o. 

INPUTS  -  Nf*>  fRI I ...  7 1  *f  « ,f,  /II l..f|.|.,-|.  NRO*Y 
URI  I...*)**.,- -t.  Ul  l...4l*0.  ,9.  ,0.  ,-9. 

OUTPUTS  -  LA*7  AAII...8l*l.,0.,0.,|,,-.)ll!),-.A,.||t| I  ,-.«*** 7 
ERR  *0. 

INPUTS  -  Nf**  fRI  l...«l*f. N.  fill ...«>•!. ,-l., •*.,«. 
NRU*f  UR  I I...B I *0., »». 4  AO., 19. 

Ull  l...*l*n.,  9. ,0. , -9. ,  0.,- I. ,  0.*  I. 

OUTPUTS  -  LA* I  ERR *0. 

AAI l..<lfl*l.,0.,0.,l.,-.)l  lb A,  .FIOiB,  .OTCbl ,- . 7107*, 

•006AA,-. 1979b,  .01  FIS,  .77*1* 


C  PROGRAM  FOIL UMS  BELUm 

c 

c 

C  Oummy  DIMENSIONS 

c 

DIMENSION  TRIf  l.f  Ilf  I,  UR  If  I.UlIf  I,  aai  f  I,  SPACE  If  I 

c 

C  BRING  IN  The  SCALAR  PARAMETERS 
C 

L  *Nf 
N*NRU 


C  SET  UP  THE  COMRINATIUNS  neeoed 

c 

"*l/N 

LA*M«| 

LN»L*N 
KN*N*N 
NNI*NN«I 
ll*l*t 
I  A  I  •  ) 

I  AM* I A|#ll 
IB  I ■  I  AR*LL 
IBR* IRI *LN 
ISP* I 8R*LN 
I  UR | • I RR 
lull  •  |R I 
I AR | ■ | AR 
I  A I  I  •  I  A  I 
C 

C  SET  UP  RIGHT  SIDE  OF  SIMULTANEOUS  EQUATIONS 
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BLANK  PAGE 


SNLPAC 


paogaan  listings 


SNLMC 


IHGf  tl 


IPACC  II 


c 

CALL  *UlA»  lUA.N.l.SPACEl  I AA  I } 

CALL  NAtAA  1 U I , N. I . S * ACE I I » I  l» 

C 

C  SET  UP  LEM  SlUt  OF  SI*UL  TANEOUS  EOUATIONS 
C 

00  10  1*1.  * 

CALL  01  ABAC  UA,2l,t,N.SPACEII0«II.SPACEIIUU.SPACHIAAt», 

I  SPACEMAN)! 

IUA 1*1  AW | 

IU 1 1 • I A 1 1 
I  AA  t  *  I  AA  I  »LN 
10  |A|  I*  I  At  I 

2000  C  I . 

C 

C  SOLVE  (HE  EQUATIONS 
C 

CALL  SIME4C  IL.L.N.SPACEI IAA I. SPACE! IAII,S®ACEI I AA I . SPACE  I IMI I . 
I  SPACE .SPACE  I ISPI.FAAI 
C 

C  IAAnSPOSE  The  OUTPUT  A SO  CHANGE  INC  SIGN  OF  OUIPUI 
C 

CALL  NAYAA  I  SPACE  I  I AA  I ,  l ,  N,  A  A  IAN  I  1 1 
CALL  HOVAEV  ILN, I.AAINNI I, I.AAINNII.-I  I 

C 

C  Put  AN  I  DENI  |tV  MAIAII  AT  IME  BE  GINNING 

C 

CALL  NOVAE V  INN.O.O.I.AA.ll 
CALI  NOVAE V  IN. 0. I . ,N* I . AA. I  I 

C 

C  THAMS  Ail 

c 

AEtUHN 

ENO 
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13.  Abstract  (Continued) 


provide  an  analytic  factorization  of  a  multi- channel  auto¬ 
correlation  in  terms  of  invertible  wavelets.  In  addition 
the  autocorrelation  may  be  approximately  factored  by  a  re¬ 
cursive  least-squares  algorithm,  or  by  a  projection  tech¬ 
nique. 

Of  the  factorization  methods  available,  the  r  *- 
cursive  algorithm  is  the  most  efficient  and  is  therefore 
extended  to  include  the  more  general  problem  of  signal 
shaping  in  the  presence  of  noise. 

Finally,  as  an  illustration,  the  problem  of  de¬ 
signing  a  finite  optimum  two-dimensional  band-pass,  band- 
reject  filter  is  solved  and  the  characteristics  of  a  few 
particular  realizations  of  such  filters  are  presented.  (U) 
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